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We study in detail static spherically symmetric solutions of non linear Pauli-Fierz theory. We 
obtain a numerical solution with a constant density source. This solution shows a recovery of the 
corresponding solution of General Relativity via the Vainshtein mechanism. This result has already 
been presented by us in a recent letter, and we give here more detailed information on it as well as on 
the procedure used to obtain it. We give new analytic insights upon this problem, in particular for 
what concerns the question of the number of solutions at infinity. We also present a weak field limit 
which allows to capture all the salient features of the numerical solution, including the Vainshtein 
crossover and the Yukawa decay. 

I. INTRODUCTION 



_^ Motivated by the wish to find consistent large distance modifications of gravity, theories of "massive gravity" have 

^^ , recently attracted some attention (see e.g. [l| for a review). Starting from the unique consistent quadratic theory for 
a Lorentz invariant massive spin 2, the Pauli-Fierz theory (henceforth PF) [^, it is easy to build a simple non linear 
completion at the price of introducing a non dynamical external metric. We will consider here such a non linear Pauli- 
^^ Fierz (henceforth NLPF) theory that we will define precisely below. One of the crucial issues to be answered in models 
of massive gravity is how to recover metrics sufficiently close to the ones obtained in Einstein's General Relativity 
(henceforth GR) to pass standard tests of the latter theory, while at the same time having significant deviations from 
GR at large distances. A first major obstacle in this way is related to the so-called van Dam-Veltman-Zakharov 
(vDVZ) discontinuity Q. Indeed, the vDVZ discontinuity states that the quadratic Pauli-Fierz theory does not have 
linearized GR as a limit when the mass of the graviton is sent to zero. In particular, the quadratic Pauli-Fierz theory 
^ ' leads to physical predictions so different from those of GR, such as those for light bending around the Sun, that one 
can rule out this theory based on the simplest tests of GR in the solar system. However, soon after the discovery of 
the vDVZ discontinuity, it was realized that the discontinuity might in fact disappear in nonlinear Pauli-Fierz theories 
tH- ' [j. This is because a careful examination of static spherically symmetric solutions of those theories by A. Vainshtein 
• ' showed that the solutions of the linearized nonlinear Pauli-Fierz theories (i.e. those of simple Pauli-Fierz theory) 
^— ^ . were only valid at distances larger than a distance scale, the Vainshtein radius Ry, which goes to infinity when the 
f— ^ ■ mass of the graviton is sent to zero. On the other hand Vainshtein showed that there existed a well behaved (as the 
mass of the graviton is sent to zero) expansion valid at distances smaller than Ry, this expansion being defined as an 
expansion around the Schwarzschild solution of GR. What Vainshtein did not show is the possibility to join together 
those two expansions as expansions of one single non singular underlying solution, as underlined in particular by 
Si^ ■ Boulware and Deser soon after Vainshtein work first appeared [5|. The vDVZ discontinuity and its possible cure a 
?H ' la Vainshtein reappeared later in more complicated (and probably more realistic) models. Those include the Dvali- 
Gabadadze-Porrati (DGP in the following) model [S] and its sequels [7], which have attracted a lot of attention in 
particular as far as application to cosmology is concerned |8l-[l3|: the "degravitation" models [lj|, the Galileon and 
k-Mouflage models [iJl, the recent models of Refs. [15l - [l7| . as well as possibly the Hofava-Lifshitz theory [18[ (see 
e.g. (H). 

The search for an everywhere non singular asymptotically flat solution featuring the Vainshtein recovery has proven 
to be a difficult task. E.g. Damour et al. found in Ref. [20| by numerical integration of the equations of motion, 
that singularities always appeared in static spherically symmetric solutions of the kind considered by Vainshtein and 
concluded that the Vainshtein mechanism was not working. This issue was recently reexamined by us, and in Ref. 
|2l| we reported briefly (disagreeing with the results of Ref. [20]) the numerical discovery of a static, spherically 
symmetric, asymptotically flat solution of massive gravity showing for the first time the Vainshtein recovery of GR. 
The aim of this work is to provide more details on this solution and on the way it has been obtained. Before doing so, 
let us however underline that we do not call for a realistic use of the kind of NLPF theory presented here (see however 
[iTl [22|y. Indeed, such theories are believed to suffer from various deadly pathologies, like the the "Boulware-Deser 
ghost" Q and the related strong coupling. Rather, we use the NLPF as a toy model to investigate the issue of the 
success or failure of the Vainshtein mechanism. 

This paper is organized as follows. In the next section |lTl mostly introductory, we first introduce the NLPF we will 
be looking at (subsection lll Ap . the ansatz considered (subsection lllBj) . the Vainshtein mechanism in brief (subsection 



IIIC|) . and then with more details, together with the so-called Decoupling Limit (DL) that played a crucial role for 
the numerical integration of the field equations (subsection III D|) . Then (section HIT]), we give some analytic insight on 
the putative solution, starting first from standard perturbation theory (subsection IIII Ap yielding in particular a large 
distance expansion, presenting later expansions valid at intermediate and small distances (subsections IIII Cj and llll Pp . 
In this section we also discuss with some details the crucial issue of the number of solutions at infinity, arguing that 
there exists infinitely many such solutions which are not captured by the standard perturbative expansion (subsection 
IIIIBI) . In the next section ITVl we introduce a new approximation able to capture both the Vainshtein recovery and the 
large distance Yukawa decay. Section |V] is devoted to the detailed presentation of our numerical results. We conclude 
in section IVll Three appendices contain more technicalities. 



II. GENERALITIES ABOUT STATIC SPHERICALLY SYMMETRIC SOLUTIONS 

A. Action and covariant equations 

We will consider the theory defined by the action 



S = / d^x.p-g ( -^R, + Lg ) + S„.t[/,5]- (1) 
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In the above action, Rg is the Ricci scalar computed with the metric g^jy, /^i/ is a non dynamical metric, Mp is a 
mass scale, and Lg denotes a generic matter Lagrangian with a minimal coupling to the metric g (and not to the 
metric /), and Sint\jT(^ is an interaction term with non derivative couplings between the two metrics. Following the 



notations of Damour et al. 23[ , these interaction terms are in the form 
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^-mHll j d\V^^\gJ) ^ -\n^Ml j d^x,/^v'^^\^-H) (2) 



with V^°'\g,f) = \J—g ^(''^(g ■"■£) a suitable "potential" density. There is much freedom in the choice of these 
potentials. In this paper we will mostly concentrate on the following term (henceforth called AGS), considered in 



particular by Arkani-Hamed et al. in Ref. [24 



S„u = -im^AfJ. J d\ V^ H^.H^r {g^'^g"^ - g^'^g'^l , (3) 

where g'^'^ denotes the inverse of the metric 5^^, and H^^, is defined by 

We will also comment on other forms for the potential. In the equation above, m is the graviton mass and Alp the 
reduced Planck mass, given by 

Mp^ = SttGat, 

in term of the Newton constant Gn- The theory described by the action ([Ij-Q is invariant under common diffeo- 
morphisms which transforms the metric as 

g^^ix) = df,x"'{x)d^x'^{x)g'„^ (a;'(a;)) , 
UAx) = d^x"'{x)d,x'^x)f'„^ {x'{x)) . 

A crucial property of the potential ([3]) is that when g is expanded to second order around the canonical Minkowski 
metric 7]^^ as 5^^ = r/^j/ + /i^j/ and / has the canonical Minkowski form ry^^ , the potential at quadratic order for /i^i^ 
takes the Pauli-Fierz form. The equations of motion, derived from action ([1]), read 

MlG^, = (T^. + Tji,) , (4) 

where G^^ denotes the Einstein tensor computed with the metric g, T^i/ is the energy momentum tensor of matter 
fields, and T^^ is the effective energy momentum tensor coming from the variation with respect to the metric g of the 
interaction term Smt- T-^^ depends non derivatively on both metrics / and g and is defined as usual as 

TjiAx) = / r uu( ^ Sint{f,g)- 

V^g "g \^) 



A simple, but non trivial, consequence of equations (jj]) is obtained by taking a g-covariant derivative V of both sides 
of the equations; one gets, using the Bianchi identities and the conservation of the matter energy momentum tensor, 

VT^, = 0, (5) 

the constraint 

WT^^ = (6) 

which the effective energy momentum tensor should obey. This equation will turn to play an important role in the 
following. 

B. Static spherically symmetric ansatz and boundary conditions 

In this paper we study static spherically symmetric solutions of the theory (Il])-(l3]) and we take the non-dynamical 
metric / to be flat, i.e. to parameterize (possibly only part of) a Minkowski space-time. Note that it is possible 
to consider / to be dynamical as well (see e.g. [20, 22, 23, 25]), for simplicity we will not consider this possibility 
here, and consider the flatness of / as a prerequisite of our model. Hence, we consider static spherically symmetric 
configurations using the following ansatz 



U.dx'-dx- = -df + ( 1 - ^^^^ ] e-^(^)di?2 + e-^<^^^R^dn^ 



Rf/{R)\ „_y(fl)jD2 , „-mW)d2jo2 *- '' 
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where i/. A, /i are unknown functions and dil^ is the metric of a unit 2-sphere 

dn^ = d9^ + sin^ edip'^ . 

This ansatz, where both metrics are diagonal in a common gauge, is not the most general one (see |25| ) however, it is 
the one relevant for the study of the Vainshtein mechanism of massive gravity which is the aim of the present work 
(see e.g. [20[). It is easy to check that for the ansatz (O, the tensors G^^, and T-^^ are both diagonal (the diagonal 
character of T^^^ follows from the fact its components are made by taking products of the diagonal matrices f , f^"*", g 
and g""*"). The tt and RR components of ([4]) and the non-trivial part of Bianchi identities ([6]) read 

^ + ]^ (1 - ^') = 8^G^ (T^« + Pe^) , (8b) 

-^^^%^« = «' ^«^) 

where the source energy momentum tensor T "^ is assumed to have the perfect fluid form 

T;" =dia.gi- p,P,P,P), 



with total mass 



The matter conservation equation reads 



fl0 

.2 



M = / 47ri?^ p dR. 



P' = -j{p + P). (9) 



In the integration of the fleld equations, we will only consider cases where the source has a constant density inside its 
radius. For the potential ([3]), we have 

Tf, = mHll, ft , T^^ = mHlj> fn, S/^Tf^^ = -mHlj>fg , 



where we use the following notations [2C 



{36"+" + e^ - 2e'') ( 1 - ^j + e^ (2e^ - e") - 3e^+^ (26^+" + e" - 2e'') 



Ir^ 



-z-' — 2/j, 



(36^+"^ - e^ - ie") (l - ^ j + e^ (2e^ + e'') - 3e^+^ (-2e^+'' + e^ + 2e'') 
X [8 (e^ - 1) (3e^+'^ - e'' - e"^) + 2R ((3e^+'^ - 2e'^) (A' + V - j/') - e^ (A' + V + v')) 



(10a) 



(10b) 



(10c) 



-i?- 



((3e^+'^ - 2e") (a'/i' - 2^" - //V + (^')') " e^ (^V ^ 2/i" + ^l'v' + (//)') - 2e'' (^i')') 



Besides the components (|8ap and (jSbp . the field equations (|4]) have non trivial 69 and </j(^ components. However, it is 
easy to check that the set of equations ([HI), (HI is equivalent to equations ^ and ([5]) once the ansatz ([7]) is chosen. 

The system ^, ([9]) is a system of ODEs that requires five boundary conditions (respectively the values of v, A, ^, 
/i' as well as P at some (finite) boundary) to be integrated along the radial coordinate R. One is also looking for an 
asymptotically flat solution such that v and A must vanish at large i?, and one should also require that the solution 
is non singular in i? = 0. To fulfill this condition and avoid a conical singularity at the origin, we will require that 



The solution found will also be such that 



A(i? = 0) = 0. 



^i\R = 0) = 0. 



(11) 



In the following subsections we first recall some properties of the so-called Vainshtein mechanism which, if valid, 
opens a way to recover solutions of General Relativity in massive gravity for small graviton mass. Then we recall 
some of our previous results on the so-called decoupling limit of the theory considered, which played a crucial role for 
the integration of the field equations. 

C. Vainshtein mechanism in brief 

The (quadratic) Pauli-Fierz theory is known to suffer from the van Dam-Veltman-Zakharov (vDVZ) discontinuity, 
i.e. the fact that when one lets the mass m of the graviton vanish, one does not recover predictions of General 
Relativity. E.g., if one adjusts the parameters (namely the Planck scale) such that the Newton constant agrees with 
the one measured by some type of Cavendish experiment, then the light bending as predicted by Pauli-Fierz theory 
(and for a vanishingly small graviton mass) will be 3/4 of the one obtained by linearizing GR ^]^. One way to 
see this is to consider solutions of equations of motion Q which are static and spherically symmetric and which 
would describe the metric around a spherically symmetric body such as the Sun. To do so, using the ansatz ([7]) 
is especially convenient because in this form, the g metric can be easily compared with the standard Schwarzschild 
solution. If one tries to find a solution expanding in the Newton constant Gn, as we recall in the next subsection, 
one finds immediately the vDVZ discontinuity appearing in the form of a different (m independent) absolute value of 
the coefficients in front of the first non trivial correction to flat space-time in gu and gRR components (neglecting the 



^ The fact it is smaller is easy to understand: the essential difference between Pauli-Fierz theory and linearized GR comes from an extra 
propagating scalar mode present in the massive theory. This mode exerts an extra attraction in the massive case compared to the 
massless case. Hence, if one wants measurements of the force exerted between non relativistic masses to agree, the coupling constant of 
the massive theory should be smaller than that of the massless theory. But light bending is blind to the scalar sector - because the light 
energy momentum tensor is traceless. Hence, provided the two theories agree on the force between non relativistic probes, the massive 
theory would predict a smaller light bending than the massless one. 



Yukawa decay by assuming the Compton wavelength of the graviton is much larger than other distances of interest). 
Indeed, those corrections are obtained at first order in GjVj i-©- by linearizing the field equations, in which case the 
NLPF and PF theories are equivalent by definition. However, as first noticed by Vainshtein [J|, the computation of 
the next order correction in the NLPF shows that the first order approximation ceases to be valid at distances to the 
source smaller than a composite scale, the Vainshtein radius defined by 

where Rs is the Schwarzschild radius of the source. This Vainshtein radius obviously diverges when one lets m go to 
zero. It is also much larger than the solar system size for a massive graviton with a Compton wavelength of the order 
of the Hubble radius and Rs taken to be the Schwarzschild radius of the Sun. Hence, one can not rule out massive 
gravity based on solar system observations and on the results of the original works on the vDVZ discontinuity [3|. 
Vainshtein also showed that an expansion around the standard Schwarzschild solution (defined as an expansion in 
the graviton mass m) can be obtained (as we recall in the next subsection). This expansion is well behaved when the 
mass of the graviton is sent to zero, opening the possibility of a recovery of GR solution at small distances R of the 
source. Indeed, the domain of validity of this second expansion was shown to be i? <C Rv- Moreover, the correction 
found by Vainshtein to the Schwarzschild solution are non analytic in the Newton constant which can explain the 
failure of the attempt to obtain an everywhere valid solution expanding in the Newton constant. 

D. The Decoupling Limit and the Vainshtein mechanism 

The Vainshtein mechanism can be enlightened taking the so-called Decoupling Limit (henceforth DL) of the theory 
considered. Indeed this limit gives the scalings of the different metric coefficients expected to hold (for an interesting 
range of distances to b e sp ecified below), should the Vainshtein mechanism be valid. The DL has been investigated by 
us in a previous work [26], and those investigations turned also to be crucial to obtain the solutions of the full system 
presented here and in Ref. ^21i] . Hence, it is the purpose of this subsection to briefly present the salient features of 
the DL and its application to study of the Vainshtein mechanism. Instead of the way it was originally introduced 
[2J|, we will define here the DL directly in the equations of motion [26]. Indeed, it can be shown that the full system 
of equations dS]) reduces to a much simpler one in the DL defined as follows 

Mp — ^ oo, 
TO ^- 0, 
A ~ constant, 
T^u/Mp ^ constant. 



where A is an energy scale defined by 



A= {m'^Mp)^ 



This scale is associated with the strongest self interaction of the model (for more details see [2J, [26|). Taking the DL, 
the equations ([5]) become 

i?//) + ^ (12a) 

(12b) 

(12c) 

where Q is the only nonlinear part left over in the DL. As seen above, the DL amounts to linearizing Gf^ and G^pp (as 
well as the sources), keeping the part of ft and fp linear in /x, linearizing fg and keeping also the piece of fg which is 
quadratic in //. As appears above, the linearization of fg does not contain any pieces depending on /i. Q is the only 
part of the pieces left over in the DL which depends on the choice of potential y^"' (provided the theory is in the 
NLPF category). For AGS potential, Q is given by 



A' 

'r 


A 
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m , 
= - — (3/i + 


v' 
R 


A 

i?2 


= rr?[i 




A 

i?2 


-;^+w 



1 //^'' ,,^4^//' 

2 -^^'^ +^ 



Q = -^M^+M/^" + ^ . (13) 



Note that we have neglected the pressure P in the above equations (rT2|) . since this is a direct consequence of taking 
the decoupHng Umit in the matter conservation equation Q- From the three equations (|12p a single equation on fi 
can be extracted, which after integration and matching of the solution to the source, reads 






(outside source), 
(inside source), 



(14a) 
(14b) 



where Rq is the radius of the source and where we assumed that the density of the source is constant. Note the r.h.s. 
of (|14bp is constant, in contrast to that of Eq. (|14ap . 

The essence of the Vainshtein mechanism, as well as the related scalings of the solutions can be understood from 
equation (|14ap . which allows to obtain /i, as well as from equations (|12ap and (|12bp which enables to get A and i' once 
/i is known. Indeed, one can guess the existence of two different regimes for the equation (|14al) , namely the linear 
and the nonlinear regimes. In the linear regime, the nonlinear part Q can first be neglected. This regime is valid for 
R ^ Ry, and the asymptotic solutions, /2oo and Aoo read 



fJ-o 



2 Rs 



3 m?R^ 



oi4 



A. 



2R^ 
3 R 



oi4 



(15) 



From (|15p. it follows that v^o — — 2Aoo, which is a manifestation of the vDVZ discontinuity, and results, e.g., in the 
mentioned difference between light bending in GR and in PF theory. In the nonlinear regime, the linear term in the 
l.h.s. of ([Ti| is neglected compared to the nonlinear Q-tcrm. In this regime (valid for R <C Rv) two different scalings 
can be identified. 

The first scaling, henceforth called Vainshtein scaling, is found assuming that the l.h.s. cancels the r.h.s. for the 
leading term in the expansion of /z. I.e., in vacuum, one assumes that the leading term /Iq is such that 



Q(/2o) ^ RslR^- 
For Q given by ((T3)) . the leading behaviours of the metric function read 

5/2' 



(16) 



^iV 



vv 



8Rs 
9R 

R^ 
R 






-,w) 



5/2' 



(17) 



T Rs 



.Rv 



5/2' 



In this regime, one recovers at dominant order the results of linearized General Relativity, namely 9y '^ —Ay. Inside 
the source, following a similar logic, one can look for a solution to the equation (|14bp JIq.v, such that the leading 
behaviour gives Q{fiQ,v) ~ Rs/R%- We find 




|mv 



{mRy + ... , 



(18a) 
(18b) 
(18c) 



where By and vq are constant of integration, to be fixed by matching the solutions outside the source, Eq. (1171) . and 
inside the source. 

A second scaling is possible, leading to the restoring of General Relativity in the non linear regime [2g. This is 
what we called the Q -scaling, whose leading behaviour corresponds to a zero mode of the nonlinear operator Q. This 



zero mode is defined such that Q(/i) ^ 0, in which case the r.h.s. term Rs/R^ compensating the nonhnear piece of 
(|14ap is subdominant. The Q-scahng solution of (|14al) reads in vacuum, 



Hq = {mRvf I Ac 



Rv 
~R 



B, 



Q 



3Ac 



■In 



R 

Rv 



^Q = -% + lAQiRvm)Hn{R/R^ 
it I 



Aq 



R 



Ac 



{RvmY 



(19a) 

(19b) 
(19c) 



From p9)) it is clearly seen that General Relativity is restored for small R. Note that the Q-solution ^T9\\ is a richer 
family compared to the Vainshtein scaling ([TT]). since it contains two constants, Aq and Bq. Inside a source, it is not 
hard to guess that the leading behaviour of the Q-scaling is also given by the Q-scaling ([T9| outside the source 



/^0,Q = {mRvf Aq 



Rv 



I'e,' 



Rs^ 
R ' 



Rs^ 
R ' 



(20) 



since the leading term of Q-scaling does not touch the r.h.s. of Eqs. (fn|) . Note, however, that the subdominant terms 
are different for the Q-scalings inside and outside the source. Note also that this scaling leads leads to a curvature 
singularity at i? = and does not obey the boundary condition ([TT1) . 

In Ref. [26|, we integrated numerically the DL system (fT2|) . We found that solutions interpolating between the 
above described asymptotic regimes indeed exist. We were able to identify the following two families of solutions. For 
a given choice of parameters Rs, m and i?©, there is one solution {flN,v, '^n,v, ^n,v} which interpolates between the 
"Vainshtein" asymptotic scalings given respectively by P^ . (|17l) and (ITS)) , as follows. 



R 


R^Rq 


RQ<t:R<t:Rv 


R->Rv 




const 

const 

const X i?^ 






VN.y 

^N,V 


^8Rs/{9R) 

-Rs/R 

Rs/R 


2i?5/(3m2i?3) 
-ARs/{m) 

2Rs/{iR) 



(21) 



Note that the Vainshtein solution is unique for fixed i?5, m and i?©, since the free constant Bv in (|18p is fixed by the 
requirement of the asymptotic flatness (|15p [26]. Apart from this Vainshtein solution, we have found a one-parameter 
family of solutions {p,N.Q, i^N,Q, ^n,q} which interpolates between the solutions given by (jlSp . and the Q-scalings (J19p 
and (|20p. such that. 



R 


R<s:Rq 


Re^R^Rv 


R:> Rv 


P'N,Q 

vn,q 


-Rs/R 
Rs/R 


m^Rl^AQ X i?-2 

-Rs/R 

Rs/R 


2i?5/(3m2i?3) 

-4i?5/(3i?) 
2i?5/(3i?) 



In contrast to the Vainshtein solution, this solution is not uniquely fixed by the asymptotic flatness. There is a 
freedom in the choice of Aq in (fT9|) . so that the family of solutions corresponds to a curve in the (Aq, i?Q)-plane of 
solutions P^ . 

Although the solutions we found in DL are solutions of a simplified version of the full system, they should be valid 
in some range of distances R for the full system. Indeed, let us estimate the validity of DL solutions (see also 26*]). 
To do this we compare the terms left over in the DL system of equations (|12p with the terms we dropped from the 
original system ([S]). First of all, keeping only linear terms in A and v on the r.h.s. of dS]) requires R » Rs outside the 
source, as in GR. The second condition comes from comparison of the linear terms containing /i with those containing 
A and v in the expansion of the r.h.s. of (j8a| and (j8bp . In the DL, the linear terms containing jjl are neglected, and 
one can see that those terms can only be neglected at distances R <C m~^ for the full system. I.e. those two first 
conditions state that the DL can at best recover linearized GR and will not capture the expected Yukawa decay of 
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the solution. The third condition arises from the neglecting of terms in the nonlinear regime. Indeed, the quadratic 
terms in /i, v and A in the r.h.s. of the expansion of (j8ap and (j8bp must be subdominant for DL to be valid for the full 
system. Similarly, the cubic terms on r.h.s. of (f8c| should be small compared to the quadratic ones, left over in DL. 
This gives nothing new for the Vainshtein scaling, namely, R 3> Rs- However, this leads to a more rigid constraint 
for the Q-scaling 

R~>Rv {mRsf^^ . (22) 

To summarize, the DL solution, having the Vainshtein asymptotic at small R and the linear asymptotic solution 
for large i?, is expected to be a good approximation for the solution of the full system in the range, 

while the DL solution, approximating the Q-scaling at small R and the linear asymptotic solution for large i?, should 
be a good approximation for the solution of the full system in the range, 

Rv{T(iRsf^'''' <^R<.m-^. 

If a not too compact source of radius Rq 3> Rs is included into the picture, one can check using the dominant 
behaviour (1181) that the Vainshtein solution even remains a good approximation of a full solution until the origin! 
This will later be confirmed by our numerical integration. On the other hand, even when a large source is included, 
the DL Q-scaling solution stops being relevant for the full system at the distance R — Ry {mRs) • 

In the following, we will look for a solution of the full system (JH), and we will use our understanding of the DL 
as a guideline. As we will report, our numerical attempts to find a solution of the full system with the Q-scaling 
asymptotic only lead to singular solutions, and accordingly, we will mostly focus on the extension of the Vainshtein 
scaling to the full system, rather than on the Q-scaling solution. 

III. INTEGRATION OF THE FIELD EQUATIONS: ANALYTIC INSIGHTS 

Before turning to present the numerical integration of the full system of equations ([8|) and ([9]) , it is convenient to 
present below some results one can obtain analytically. Those will mostly be given in the form of series expansions. 
We have no proof of the convergence of those series, and we in fact believe (and will give some arguments below) that 
those series are divergent, being then asymptotic expansions. However, we found those series useful for three main 
reasons: First they allowed us to set boundary or "initial" conditions for the numerical integration with the required 
precision; second, they were also useful to check that the numerical solutions we found was not singular at the origin 
of the radial coordinate system; and third, because they help investigating the uniqueness of the solutions at infinity^. 

Let us first outline the organization of this section as well as its main results, allowing a reader not interested in 
technical details, to skip all but this introduction. The first subsection IIII Al explains with some details the solution 
of standard perturbation theory. This perturbation theory is defined as an expansion in a small parameter e where 
at each order, all the metric function A, i^ and /z are assumed to be proportional to e raised to the same power. The 
parameter e is expected to be proportional to the mass of the source (or to its density) . Solutions of the perturbation 
theory out of the source are obtained at each order in e in the form of asymptotic series both for small distances and 
for large distances. Those series are expressed in terms of the dimensionless quantity z — Rm. Hence, we obtain a 
double series expansion (into e and z). The "small distance" expansion, which is an expansion valid at small z (at a 
fixed order e), is similar to the regime ()15p of the DL, and is found (by comparing successive orders in e) to be no 
longer valid at distances smaller than the Vainshtein radius in agreement with Vainshtein's computation. This "small 
distance" series expansion depends on two unknown integration constant which have to be fixed by matching to the 
source. On the other hand, the "large distance" (large z) expansion is uniquely fixed and has no free parameter. 
However, we argue in subsection IIII Bl that this is not enough to define a unique solution at infinity, namely because 
we show that a non perturbative correction can always be added to the unique solution of perturbation theory we 



^ Similar series expansions were used in Ref. [2011 to argue that the system we investigate here has a unique asymptotically fiat solution. 
It was then argued on the basis of this, that it was "a priori" unlikely that one could match this unique solution to a solution obtained 
d la Vainshtein, expanding around the Schwarzschild solution, as will be explained in section IIII CI As we will see our results do not 
support this conclusion on the uniqueness of the asymptotically flat solution, hence removing the argument given in Ref. ^U] to explain 
the singularities found by the numerical integration done in this last reference. 



identified. This is argued to be related to the fact the series presented are not convergent series, in analogy with 
what we have shown for the decoupling limit |26| . Then, we present expansions valid inside the Vainshtein radius first 
outside the source (jIII C[) then inside it (|IIIDp . The expansion inside the Vainshtein radius (jIIIC|) . but outside of the 
source, is obtained expanding the functions around GR, in a small parameter e' which is by definition proportional 
to m^. At each order in e', one can then expand the functions in powers of R/Rs {Rs being the Schwarzschild radius 
of the source) . One thus obtain again a double expansion. This expansion, allows to compute the first non trivial 
corrections to the GR solution. Those corrections will then later be compared to the ones found by the numerics. 
One key property of the Vainshtein mechanism is related to the existence of those double expansions, to the fact the 
associated series appear not to be convergent allowing the addition of "non perturbative corrections" such as the ones 
mentioned above, and to the fact, as discovered by Vainshtein, that the solution in the "below Ry" regime features 
corrections to GR which are non analytic in the e parameter of the standard perturbation theory. Last (jIIIDj) . we 
obtain a new expansion valid inside the source, that will be later compared to the numerical solution, and also used 
by us to study the regularity of the solution at the origin. 

A. Standard perturbation theory in vacuum 

Looking for a solution of the system ([5]) , © far from a source, a first obvious way is to look for a scheme where the 
functions A, v and ^ can then be expanded as 

A = A0 + A1 + ... 

iy = VQ + vi + ... (23) 

^ = /.to +A*i + ■•• 

where A^, Vi, fii are assumed to be proportional to e'+^, e being a small parameter, and such that A^+i ^ A^, I'i+i ^ Vi 
and /ii+i ^ /ii, at least for some range of distances. This is what we will call here standard perturbation theory. 
Before going further, let us first introduce a new dimensionless radial coordinate z defined by z = Rm. In this unit, 
R = m~^ corresponds to z = 1, and R = Ry translates to z = a, where a is defined as 

a = Rym. (24) 

Using this variable, the lowest order Aq, i'q and /io are obtained solving the system (I8all8c[) linearized which reads 

A' A 1 

— + ^ + - (A + 3/1 + z/i') = 0, 

z z^ 2 



1^' A 1 

7^ ^ ^ 2 



---o(^ + 2m) = 0, 



, 2A 

z 



The asymptotically flat solution of the above linear system can be found analytically, it is given by [20|, |27| 

4b 

1^0 = -C -e-^ (25a) 

Ao = C y A + i") e-^ (25b) 

^ 26 / 1 1 \ , , , 

Mo=C— 1 + - + - e-^ (25c) 



where 

b = a^ ^ {Rymf = Rsm, 

and C is a dimensionless integration constant. The value of the constant of integration C can in principle be determined 
by integrating the system inwards and matching to the source, as can be done in a similar case in GR. There, indeed, 
an integration constant with the same status arises, when one integrate the linearized Einstein equations in the vacuum 
around a static spherically symmetric source. This constant of integration is then seen to give the Schwarzschild mass, 
once matched to the source. Here, an alternative way to proceed, is to use our knowledge of the Decoupling Limit. 
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Indeed, requiring the solution to be very close to the Decoupling Limit solution dTSl) in the range Ry <Si R -^i m~^, 
i.e. a <^ z <^ 1, fixes the constant to be 

C = 1. (26) 

In Section IV Bl we will show that this choice of constant of integration is also in agreement with the numerical 
integration. Accordingly, from here-on and until the end of this article, we will assume (except at some scattered 
places) C = 1. Note however, that when dealing with the equations in vacuum, the value of this constant is irrelevant. 

Having fixed C one notices that the solution (f25|) features the expected behaviour for the gravitational potential 
of a massive graviton: it has a Yukawa decay and the "Newtonian" part i^q is 4/3 larger than the GR result (in the 
range of distances R ^ m~^) while the radial part Ao is 2/3 smaller, a structure related to the vDVZ discontinuity. 
Our expansion parameter e, can then be defined as e = C x b, note that it is indeed much smaller than one for usual 
astrophysical sources and a graviton of cosmological size Compton length. 

To go beyond the lowest order (P5)) . the first step consist in re- writing the system of equations ([5]), separating 
explicitly the linear terms from the nonlinear ones: 

A' A 1 

l-^ + 7:(A + 3^ + z/i') = ft,^2 - Gtt,^2 , 

z z^ 2 

z z^ 2 / .^ ,^ 



z^ 2z z 



-2 9^ ^■lg.>'2- 



In the above system, we have used the expression (1101) for the functions fi and the /i,^2, and Gii^^2 stand for the 
nonlinear part of the functions fi and of the components of the Einstein's tensor in spherical coordinates. At each 
order n, the functions A„, Vn and ^„ (with n > 1) can be obtained by solving the above system where the left hand 
side consists in a linear operator acting on the unknown A„, Vn and /i„, while the right hand side is obtained keeping 
the order e"+i and consists in non linear expressions depending on the Xi^Vi and ^i (with i < n). As will be shown 
in detail in the appendix [^ the functions X-mV-a and /j,„ (with n > 1) can be expressed as asymptotic (and formal) 
expansions in large z (large distance) or small z (small distance) regimes. We now discuss the main properties of 
those expansions, starting, for reasons we will explain, by the small distance expansion. 

1. Small distance expansion 
An expansion, asymptotic at small z, can be obtained in the form 

n-f 1 oo n 

. n — k 



/^„ = z^fe^-j ^ ^ ^„^j,fe z* (z^ log 

^ ^ ^ 2=0 fe = 

/ _^ \ n+l oo n 

A„ = z4(e^) ^^A„,,,fe z^(z5 log z)"-\ (28) 

oo n 
^^'^n.i,k Z' {z^ log z) " , 



;/„ = z"* I £ 5 



' i=0 fe=0 

— ^ \ n+1 oo n 



i=0 fe=0 



Note that the form of this expansion was sketched in Ref. [2fl|, with however some important differences: the 
corresponding expansion given there is stopping at a finite value of i at each order n and did not contain any 
logarithms. As shown in the appendix |X1 it is however not possible to obtain such a simplified expansion. As detailed 
in the same appendix, two important properties of the expansion (j28p are as follows. First, when solving for the 
order n, the coefficients fJ-n.i,kT Xn.i,k and i'n.i,k are determined only after choosing two arbitrary constants. This is 
because the linear operator appearing on the left hand side of equation (j27p has two independent zero modes, (among 
which the zero order solution (1251) ) which can be expanded formally in the same way as above'^. One way to fix those 
constants would be to enforce at small z the hierarchy A^+i <^ Xi, Vi+i <S^ Vi and /i^+i ^ /i^. We have in fact checked 



The presence of logarithms in this expansion is also related to those zero modes. 
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that it is possible to do so at order n — 1. However, it is hard to go further because it is far from obvious that the 
expansions considered are convergent. The question of convergence is actually quite hard to address, since it requires 
to know the behaviour of the coefficients Xn,i.k, i^n.i,k and ^n.i,k- This can be achieved quite easily for a subset of 
those coefficients. Indeed, let us, at each order n retain the dominant term in the expansion (j28p at small z (i.e. for 
R ^ m^^). We obtain, for z ^ 1 the following expansion for /i 

M-E^"^^gfH^- (29) 

Equivalently, this is what is obtained in the limit tti — > 0. In fact, and this is the second important property we would 
like to mention here, as shown in appendix lAl the expansion ([29l) is exactly the one obtained in the DL (c/. Ref. [261], 
Eq. (4.34)). The coefficient /^n,o,n can be easily computed in the DL, and it has been shown in [26] that this series 
diverges for any z, providing only an asymptotic expansion of the solution. One can even go further: using the exact 
values of the coefficients fin,o,n obtained in the DL, one can show that the subseries of ([28| given by 

^^„+i^-(„+i).^^ (30) 

diverges for any z, despite the presence of the exponential damping e~("+i)^. Of course, the fact that this sub-series 
of ((28|) diverges does not lead to any conclusion by itself besides the fact that the expansion (|28|) cannot be absolutely 
convergent (if it were the case, any subseries would converge and one could rearrange at will the order of summation). 
In our case, there is no clear prescription for the order of summation, and the divergence of the sub-series pop does 
not mean that the full series diverges. Still, we see this divergence as a clue that the convergence of the full series is far 
from guaranteed, and that any conclusion based on these expansions should be taken with great care. For instance, 
one may question the approach followed in [20|, Sec. IV. B, which consists in building the solution of massive gravity 
far from the source and arguing on its uniqueness, starting from an expansion of the form (j28p . which is implicitly 
assumed to be convergent. Moreover, as also happens in the DL, even if the series expansion ([28| turns out to be 
uniquely defined (once the constants we mentioned have been fixed e.g. by the procedure outlined above), this does 
not mean that it defines a unique vacuum solution because it could be, if the series is not convergent, that a family 
of solutions all have the same asymptotic series expansion. We will in fact come back later to this important issue. 

Though, let us state that from a numerical point of view expansions ([25]) and (|29l) are very useful. Indeed, if one 
keeps only the relevant first orders, these expansions appear to be excellent approximations of the solutions and can 
be used to set properly boundary or initial conditions in the numerical integration. 

Moreover, by comparing two successive order in e, it is easily seen that for R ^ m~^, the expansion ((28)) ceases to 
be valid when R becomes smaller than the Vainshtein radius (i.e. whenever ez~^ becomes of order one). 

2. Large distance expansion 

In a similar way as above, an asymptotic expansion at large z can be found at each order n in the e expansion. It 
is derived in appendix |X1 and reads 







i=0 


Mn 


= e"+ie- 


/ , tJ-n,iZ , 

i — — OQ 

i=0 


An 


= e"+ie- 


in+l)z \ ' X -7* 
/ , ■^n,iZ , 

j— — oo 

i=0 


l^n 


= e"+ie- 


("+^^^ > ; ^„,.^^ 



(31) 



One can verify easily that the zeroth order exact solutions (I25p can be expressed as above. In contrast to the expansion 
(1^51) , there is here no need for logarithms nor to fix "constants of integration" . As such, this could indicate in a clearer 
way than in the previous case that, as stated in Ref. [20], perturbation theory defines a unique solution. This 
is however again not obvious because this relies crucially on the fact the above series converges. In fact, we have 
computed numerically the first coefficients of the first non trivial order n — 1 and check that the corresponding series 
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seems to have a vanishing radius of convergence^. In the next section we try to address the issue of the uniqueness of 
the solution at infinity in a different way. This analysis again indicates, in agreement with our numerical investigations, 
that the solution at infinity is not unique. 

B. Investigations of the number of solutions at infinity 

As underlined in the previous subsection, it is hard to conclude anything on the uniqueness of the solution at infinity 
from the existence of series expansions, which can be divergent. Here we would like to provide arguments showing 
that there are in fact infinitely many asymptotically flat solutions at infinity which have the same dominant behaviour 
at large distance once a constant, corresponding to the "mass" of the source in the dominant terms has been chosen. 
This non-uniqueness appears in the form of a non perturbative correction that can be added at will to the unique 
solution of the perturbation theory.^ This is in striking contrast with the General Relativity case. Indeed, in GR, 
once the asymptotic behaviour has been chosen through® v = —Rs/R and A = Rs/R, there is a unique solution 
which follows this asymptotic behaviour at infinity: the Schwarzschild solution which reads v = vgr = lii(l — Rs/R) 
and A = Xgr = —vgr- In NLPF Gravity — at least for potentials of the type (jS]) — we shall see that the asymptotic 
behaviour is not enough to determine uniquely the solution. 

To prove the existence of these infinitely many solutions, let's first assume that there exists a solution of the full 
nonlinear system of equations (|S]), a = {A, v, /2}, such that 

(T — > (To for z — > oo. 

where ctq = {Aq, I'd, /^o} as defined in Eq. (j25|) . and we imply we have fixed all the necessary integration constants such 
that a and a^ are fully specified. We then want to show there is an infinite number of solutions in vacuum having the 
same asymptotic at the infinity. We will look for these solutions a as small perturbations of the known solution a: 

a — a + Sa, 

where da = {6X, dv, Sfj,} and we will further assume that 

|<5a|«|fT|, |<5a'| « |a'|, |<5a"| « |a"|. (32) 

Taking into account (j32[) . we linearize the system of equations ([5]), for z — >■ oo, around the solution a to obtain the 
following system of equations on (5A, Sv and Sfi: 

— + ^ + liS^ + 3^M + ^V) = ^^'^, (33a) 

z z^ 2 dz 

— -^-^(5^ + 2<5m)-0, (33b) 

^ - ^ = ^ f <5^" - V - 2^^ + ^ + 25.) . (33c) 

2z z^ 3z^ \ z J 

Note that we substituted CTq instead of a in p3p . since ct — > cto at z — >■ oo and we are interested in solutions at infinity. 
Until now, we have not made any assumptions about the form of 6a, except (I32[) . If in addition one assumes that 
the r.h.s. of (|55)) can be neglected, then the solution for Sa will be given by (P5|) provided that the integration constant 
b is replaced by Sb, such that \Sb\ <C \b\. However in this case one should impose Sb = since we require ct — J- ctq at 
z — > oo, thus we do not obtain any new solution (this is, mutatis mutandis, what happens in GR). 
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One finds numerically that |/ii i//Ji i+i| oc |i| indicating a vanisiiing radius of absolute convergence 

^ Note that we found similar result in the DL case [26] , as recalled in Sec. IIIDI Indeed, in the DL case, these infinitely many solutions 
allowed to accommodate several types of solutions at small distances: the Vainshtein solution I I17I I and the Q-scaling solutions II19I I 
which depend on a free constant Aq. These solutions also appear as "non perturbative". It is interesting to note that we can arrive 
to the same conclusion (that there is an infinite number of solutions at the infinity with the same asymptotic behaviour) if we consider 
the weak-field equation approximation (see below, Sec. lIVI l as a starting point for this proof. The general structure of the additional 
solutions remains the same as for NLPF (see below in (|4TJ), with, however, different exponent and pre-exponent factors. This difference 
is due to the fact that when getting an analogue of II33I I from the weak-field equations I I59I I. some terms will be missing compared to 
ll33l l. 

® We recall that in GR, the vacuum solution depends on one integration constant, which can be identified with the Schwarzschild radius 



13 

Instead, let us try to find another solution such that (some) terms containing e~^ are not negligible in (p3| . In 
particular, we would like to keep the term containing S^" . Introducing a new variable C as follows, 

the system of equations ([55]) can be rewritten as, 

^'' ' '' +Us>^ + ^S, + Clo,CS,)=J^, (34a) 



logC (logC)2 2V . -, .. ^j 3^^^g^ 

^^^ ^^ -l-iSiy + 2St,) = 0, (34b) 



logC (logC)2 2 

^ „ J^ ^ ^^ fc^j^ _ 25^ + -^ + 2<5.^ , (34c) 

2iogc vc 3ciog2c \^ logc ; 

where dot denotes the derivative with respect to (^. In order to simplify the system of equations (1341) , let us make 
further assumptions about the functions 6\, Sv and Sji at (^ ^ co: 

\6ly\<$:\s^i\<$:C\Sfl\<$:C^\S^l\<$:C^\6^i\, iogC|<5A| < CI-^AI < CI-^Aj. (35) 

The equations (|M)) then take the simpler form 

5A + i (log 0^(5^ = 0, (36a) 

6.- -^-'^6,^0, (36b) 

^'^ + ^ {^^^ - Clog C S'ly) = 0, (36c) 

These equations can be combined to give one differential equation on /z: 

CSii + 5ii - ^(logO^^M - ^'^M = 0, 
which — applying the assumptions psp — can be rewritten as 

5^i~^^{\ogCf5^x = o. (37) 

The asymptotic form of the decaying solution of Eq. ([57]) can be easily guessed: 

5/. = Foo(C)expU3y^logCj, for C -> «^, (38) 

where -Foo(C) is a slowly varying function compared to exp I — it^^^ ] that cannot be fixed through this leading 
behaviour analysis (c/. Appendix [Cj for a more detailed approach based on series expansions). Note that a companion 
growing mode can also be considered 

5/. = ^i,s™-)(C)expf3y|logCJ, for C ^ oo, (39) 

which is fixed to zero by the conditions p2p : this growing mode can play an important role while integrating the 
equations of motion a la Runge-Kutta, since it can blow up very rapidly (as exponent of exponent) if sourced 
by numerical errors. This explains the difficulty one encounters while trying to integrate the equations of motion 
for distances beyond Ry for the potential ^. In the following, we will discard such a term since it contradicts 
the boundary conditions at infinity, but one should keep in mind its possible existence when solving the equations 
numerically. 
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Prom ([36)) one can find the other functions at C ^- oo, 

S\ = _Foo(C)^i^|^exp Us^llogCJ , (40a) 

6^ = -F^(C)^i=exp Usy^logCj . (40b) 

Finally, one can check that the assumptions (|5^ and p5l) are satisfied for the solutions l\'38\i. (HO]) as C -> oo (z — > oo). 
In terms of the variable z, equations ([55]) and (jlHl) read 

(5// = i^oo(-2)exp f ^ze^/^j, (41a) 

a = -F^(z)yexp('--|ze^/2V (41b) 

5^ = -Fo,(z)^cxp (^-| - 7^ "^'^') ■ (41c) 

Note that it is possible to complete this leading behaviour analysis through a series expansion of the solution; we refer 
the reader to Appendix [Cj for more details. 

Thus we have shown that apart from the solution A, D and /2, there is at least one family of asymptotic solutions 
at z ^>- oo, parametrized by the function Fq^, and given by 

X^\ + 6X, 
u = D + Sv, 
u = D + Siy. 

Note, that the difference between solutions decays extremely fast — as a double exponent — as z — )■ oo. 

C. Expansion inside the Vainshtein radius and outside the source 

Let us now discuss what happens inside the Vainshtein radius, but outside the star. This will prove useful for 
different reasons. Wc will indeed recall how one can obtain an expansion around the GR solution, and hence compute 
the corrections to GR first found by Vainshtein. This will later be compared in section FVBI to the numerical solution. 
We will also comment on the role of the integration constants appearing in the expansion procedure. 

Let us first rely on our investigation of the DL, which provides interesting insights into the different types of 
solutions in the region Rq ^ R -^ Rv ■ It has been shown in [26] (as also recalled in Sec. IIID[) that, in the DL, two 
types of solutions exist in the region Rq ^ R ^ Rv'- the so-called Vainshtein solution of Eq. p7)) and the Q-scaling 
solutions p9)) . We are interested in obtaining, in the full theory, a solution which will be close to one of these DL 
solutions, for a range of distances as broad as possible. It has been recalled in Sec. Ill Dl that the Q-scaling solutions 
stops being relevant for the full nonlinear theory at the radius R^ ~ Ry (mRs) ■ Below this distance, the DL is no 
longer useful for understanding the theory: new nonlinear terms enter into the game and one has to rely on numerics 
to investigate the existence of solutions. We will come back to numerics in Sec. |Vl but we can already announce that 
our investigation has been negative in this direction: it has not been possible to us to find a regular solution of the 
full system which would exhibit a Q-scaling behaviour in the range R^ = Ry (mRs) <^ R <^ Ry- all our numerical 
Q-scalings solutions were developing strong instabilities around the problematic radius i?*, and integration was not 
possible for smaller distances. For this reason, we will mostly focus in the rest of this article on the Vainshtein solution, 
and try to generalize the DL version of this solution to the full system, under the form of series expansions. A first 
resolution method consists in looking for an expansion whose zero-th order is just the DL. A convenient parameter 
for such an expansion is the parameter a = Rym introduced in Eq. (j24p . Indeed, one can show [26] that after a 
rescaling of the functions and equations, the limit a — >■ in the system ^ leads to the DL equations (|12l) . The full 
solution expansion reads 

oo oo oo 

fi = a^J2ni^\R/Rv)a^'', A = a^ ^ A(f'(-R/-Rv) «'", J/ = a^ ^ i^i")(i?/i?y) a"", (42) 

n—0 n—0 n—0 
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where the zero-th order functions a^fii^'{R/Rv),a'^\Q{R/Rv),0'^VoiR/Rv) are the solutions of the DL system, 
whose expansions are given in Eq. ([TT]) . Even if the expansion P^ has the advantage to Unk exphcitly the solution of 
the full system to the DL, it does not allow for an easy comparison with the GR solution. For this purpose, another 
expansion is possible as was originally proposed by Vainshtein [4|. 

Indeed, it is possible to look for a solution as an expansion in powers of the (square of the) graviton mass m^ dealing 
with the mass term as a perturbation around massless GR. Let us define this expansion as 



fJ- 



E 

fc=0 



fi^^''\R/Rs) m 



2fe 



X^y\'^^\R/Rs)m'\ V 



k=Q 



E(m) 
''k 
fe=0 



{R/Rs) m 



2k 



(43) 



By definition, it is such that the lowest order is equal to the Schwarzschild solution 

2 



A^™) = -i.^") = - In (1 



R 



R 



El 
R 



(44) 



Following Vainshtein, let us now briefly sketch how to obtain the first other non trivial terms of the expansion. The 
function ^q can be found using the Bianchi identity (j8c|) . where A and v have been replaced by their GR expression 



\(™) ,,('"). 



R-^' 
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= 0. 



(45) 



One can look for a solution for /in as an expansion in powers of {Rs/R), valid far from the source, where R 3> Rs- 



Ai6 



In this regime, the functions Aq 



(m) (m) 



^0 



go to zero, meaning that it makes sense to linearize the Bianchi equation in 



Aq and j^q in order to find the dominant behaviour of ^j-q ■ In addition, if fij"" is assumed to also tend towards 

zero, it is consistent to only keep the terms of lowest order in this function. Since there is no linear term in /ip i^ 
the Bianchi equation, this lowest order is nothing else than the quadratic term Q defined in Eq. (I13p . The Bianchi 
equation (j45p then reads, at dominant order 



i?2 



2R 



Q(Aio) 



Using the asymptotic behaviour far from the source of the Schwarzschild solution (1441) . i.e. the Newtonian regime, 
one gets the equation 



which is identical to the DL equation ([T 
in R, one finds 



If one assumes that the function /ig *^an be expanded as a power series 
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(46) 



The higher powers in R can then be found order by order in R [20|. Let us now compute the first order terms. The 

functions A^"* and v^ are computed through the Einstein's equations (|8a|) and (|8bl) . In the limit Rs <^ R <^ m^^, 
these equations read 
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Using the behaviour (|46p , the leading behaviour of the first order functions can be found to be 
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(47) 



while the function /ij^™ is obtained using the Bianchi identity and reads 
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(48) 



Comparing the zero-th order terms, given by Eq. pij) and Eq. P5|) . and the first order ones, given by Eq. P7)) and 

Eq. (|48|) . we reach the conclusion that the expansion in powers of m? is only valid for R < Ry = {mT^Rs) ■ Hence, 
the functions ([47|) give the first non trivial corrections to the Schwarzschild solution for distances below the Vainshtein 
radius. In section |Vl it will be seen they agree very well with the numerical results. Notice also that those corrections 
are non analytic in the expansion parameter e of the perturbation theory of section IIII Al It is also interesting to 
notice that the expansion presented above and the "small distance" perturbative expansion of subsection IIII A ij both 
yield the same "Vainshtein radius" Ry as a limiting distance of their respective domain of validity. 

One may wonder which of the expansions (|42l) and p3)) is the most appropriate for describing the solution of the 
full NLPF theory. There is no simple answer to this question since, eventually, all the orders are needed to fully 
recover the solution. Still, a simple but instructive calculation can be done: one can check that for R < (Rym)^ Ry , 
the correction to the Newtonian potential v^ ^ —Rs/R due to the GR nonlinearities {i.e the first correction to 
vn in Vq ) dominates over the correction due to the DL nonlinearities (i.e the first correction to v^q in a^v^ ). 
In mathematical terms, this reads (Rs/R)"^ ^ (mR) ^jRs/R- In particular, this means that for sources of radius 
Rq <^ {Rym)^^'^Ry (which is usually the case for sources we considered in our numerical investigations), the GR 
solution is a better approximation inside the source than the DL solution. We were able to confirm this simple analysis 
numerically (c/. Fig. [3]). 

Note that here, for the sake of simplicity, we have followed |20| and discarded any integration constant appearing 
in the process of obtaining the successive terms (beyond the GR part (|44p ) in the expansion (j43|) . Such constants 
appear in particular at each order via a double integration of the Bianchi identity. If there is however a solution with 
such an expansion that can extend from i? = to i? = oo, there is no reasons for such constants to vanish, and they 
should be kept explicit until one does the matching with the source. 

D. Expansion inside the source 

Let's now turn to the interior of the source, that we assume to be of constant density p and of radius Rq. In order 
to gain some understanding of the solution near the origin, it is possible to expand the solution as 

oo oo oo oo 

n— n— 1 n^O n—0 

where we have defined the rescaled variable 

R 

Note that the expansion of the function A starts at n = 1 (in agreement with condition ([TT|) ). These series depend on 
3 integration constants, that can be chosen to be fx^' = p,{R = 0), v^' = v{R = 0) and P^' = P{R = 0). The other 

coefficients A'So' '^n>oi ^n>o ^^^ then functions of these integration constants /iQ , i'q ,Pq and the mass ra of the 
graviton. These coefficients can be found solving order by order in x^\ a few of them are given for the AGS potential 
in Appendix [Bl 

IV. WEAK-FIELD APPROXIMATION 

As we have just seen, different expansion schemes can be used in order to get understanding of the putative solution. 
However, those schemes all fail to describe correctly all the key features of the solution. More specifically, the DL 
scheme (section [II Dp describes well the Vainshtein crossover, but misses the Yukawa decay. The perturbation theory 
approach f section [HI Ap catches the Yukawa decay, but misses the Vainshtein crossover, while the expansion around 
GR misses both the Yukawa decay and the Vainshtein crossover (section IlIICp . The purpose of this section is to 
introduce an expansion scheme, called here weak-field approximation, that will retain both the Yukawa decay and the 
Vainshtein crossover. 
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We expand the full system of Einstein equations together with Bianchi identity in powers of A, v, fi, and assume 
the smallness of each function fj,, v and A as well as of their derivatives: 

{A, V, /i} < 1, {i?A', Rv' , Rn'} < 1, R^fi" < 1. (50) 

Expanding ([5]) and using our assumptions (j50p we obtain, 

^ + :^ = -^ (^ + ^^ + ^^' + ^("')) + w ^^'"^ 

^-4 = ^(- + 2m + 0{X') + 0{R'^P)) , (51b) 



X-^ = ~R^(^ + ^ + ^]+ 0{v^) + 0{\'^R^) + 0{v>!R) 



R;/ _ 

+ 0({i^,A'}x {^,/^'i?,^"i?2}). (51c) 

Note that P in the r.h.s. of (|51bl) disappears as a consequence of the conservation equation © and the assumption 
of the weak field regime, similar to the DL case. One can recognize, in the r.h.s. of Eq. (|51cp . the non- linear DL part 
responsible for the Vainshtein mechanism. Note that so far we did not make any assumptions about relation between 
A, V, fi or between functions and their derivatives, and this is the reason why we kept the quadratic terms in fi on 
r.h.s. of (l5Tc)) . 

To further simplify the system of equations (j5ip we assume that the quadratic terms ^ i^^ are small in comparison 
to the linear terms in r.h.s. of (|51ap . i.e., 

v^ <S^ max{A,/i} , 

and similarly that the term ~ A^ in the r.h.s. of (I51b[) is negligible compared with i/. We also assume that the 
quadratic terms containing /i and its derivatives are negligible compared to the linear terms containing fj,. It is 
important to note that the last assumptions eliminates the quadratic terms in Eqs. (15 lap and ()51bp . while we have to 
keep the quadratic terms in /i in (I51cp . since the linear term in /i is absent in this equation. Note, that if we assumed 
at this step that ^ ^ v ^ X, then the quadratic terms in {/i, /i',/i"} in the r.h.s. of (|51cp would be dropped. It is 
important, however, that we do not assume any relations between {/x, ^' , fj,"} and {ly, A, v' , A'}, and thus we keep these 
terms. This will be the key feature of the expansion scheme introduced here. As a result we obtain the system 



A' A m 

R^R^ 

v' 



2 

(A + S^i + R^JL) + SttGatp, (52a) 



A m , ^ , ^ 

^-;^ = -0(m,m',m"), (52c) 

with Q given by ([T3|. It is then convenient to use once again the rescaled radial coordinate z — Rm introduced in 
subsection IIII A\ such that we obtain from (f52|) . 

A' A 1 

— + ^ = --{X + 3fi + £,!./) +p, (53a) 

z tJ 2 

""' ^ ^{v + 2^i), (53b) 



z z2 2 
A __i/ 

"~ Yz 



2 ,- = Q(/i,A*',A^"), (53c) 



where 



and we introduced p = SttGnp/tit?- This system of equation corresponds to the one left over in the limit considered 
here. As will be seen, it has the required properties. 
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First, note that it is possible to rewrite this system as a single ODE on /i and two algebraic relations between A, v 
and [i and its derivatives. Indeed, expressing A from (|53bp in terms of v, v' and fi and substituting it to Eqs. (j53ap 
and (|53c|) we find two equations containing only v and /i (and their derivatives): 

u" + — --{z^ + Q)v- -zfj.' --{z^ + 3)fi = p, (55a) 

-^+'^+^L + Q = 0. (55b) 

Taking the derivative of (j55bp with respect to z and substituting the resolved i^" into (|55aP we obtain a system of two 
equations: 

^^' - i (z^ + 6) ^ - 1 (z^ + 3) /. + ^ (4Q' + 3m') - p, (56a) 

-|; + ^^ + /i + 0-0. (56b) 

From the above system of equations (1561) we can find i/' in terms of /i and its derivatives and v in terms of fi and its 
derivatives, 

i/(z) = J", 
iy'{z) ^ zT+2ziQ + p). 



where we denoted 



^^- 3(^2 + 2) [^ ^^^' + 3^') + (^' + 3) (4^ + 3/') - 2p] (58) 



Thus we arrive to the following single ODE on fj,, 



1 -T^ 

-— ^zF+2z{Q + p). 
dz 



The above equation can be rewritten in a closed form. 



(2 + z^) (4Q" + 3/i") + - (4 + z^) (4Q' + 3/i') - (z^ + 4) (4g + 3^) - Q {z^ + 2)' 

= (2 + z2)^_2(4 + z2)p. 



(59) 



Note that Eq. ((59| is of the fourth order, meaning that the number of boundary conditions needed to solve the 
equation is four. This is in agreement with number of boundary conditions needed to specify for Eqs. ()52p . Eq. (|59p 
can further be simplified under appropriate assumptions. First of all, let us recover the solution (f25|) for the linearized 
equations. To do so, we drop in (j59p all non-linear terms, i.e. all the terms containing Q and its derivatives, as well 
as the source terms (containing p and p') to get 

(2 + z^) //' + ^ (4 + z^) p' - (z^ + 4) M = 0. 

The non-growing at the infinity solution of the above equation reads 

_ Ce-^(z2 + z + l) 

/i ;^ , 

z-' 

which coincides with Eq. (PSJ) and catches the Yukawa decay. 

Another limit can be obtained assuming in ((5^ as follows. We assume Qz <C Q' and pz ^ /i' together with z <C 1. 
Neglecting sub-leading terms in ((59)) . one obtains, 

(4g" + 3/i") + - (4Q' + ip') = ^ - 4p. 

z z 

The above equation can be integrated twice to give, 

2Q+^p=^(co+J z^pdz\ + Ci . (60) 



where in the last expression we also assumed, 
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z / pzdz \ dz <^ z pdz 

The constant of integration Ci in (pO]) must be set to to get correct behaviour of p at infinity (otherwise the solution 
is non-flat at large distances). Similarly, Cq = 0, by matching p inside the Vainshtein radius and outside the source. 
Restoring physical units, R = z/m, one can see that we arrived to the previously found DL equation (|14p catching 
the Vainshtein crossover. Furthermore, the analytic solution of the weak field system (|53p will be seen to reproduce 
very well the numerical integration (in the expected range) in the following section. 

V. INTEGRATION OF THE FULL SYSTEM: NUMERICAL RESULTS 

We have integrated the equations of motion ^ and (JH) using two different approaches. The first one is based on a 
resolution by relaxation, which turns out to be very well adapted to the problem. This method and the related results 
are presented respectively in the subsections IV Al and IV Bl The second approach, which is described in the subsection 
IV CI consist in a shooting method a la Runge-Kutta. Both methods perfectly agree, demonstrating the robustness of 
our numerical investigations. A reader mainly interested in the numerical results can read only subsection IV Bl which 
is self-contained. 

A. Relaxation approach: method and boundary conditions 

The basic idea of the relaxation method consists in choosing an initial configuration for the solution, and then 
deforming it until the equation we want to solve is satisfied with a good enough precision. The initial guess is of 
prime importance, since it will determine whether the iteration converges or not. In the case of the NLPF equations 
of motion ([5]) and ^, a clever guess is needed: indeed, starting from a naive guess {e.g all the functions to zero) does 
not lead to a converging iteration. Here is where the DL plays a crucial role in finding a solution of the full system: 
we can use the DL solution as a guess for the full solution. Note that the DL being simpler than the full system, the 
DL solution itself is not too hard to find (no specific guess is needed) . 

The advantage of the relaxation approach is that boundary conditions can be easily enforced. In particular, this 
method allows to fix boundary conditions at different points, for example at the two ends of the interval on which we 
want to solve the equations. This flexibility in the boundary conditions fixing makes this method very well fitted for 
static solution studies in which we want to impose conditions both in the center and very far from the source. 

In our case, we worked with the dimensionless variable 

C = R/Rv, 
as well as with the rescaled functions w, u, and v defined by 



w 


— 2 

= a /i. 


V 


—4 

— a i^, 


u 


= a'^X. 



We fixed the boundary condition at ^ = R/Ry — and at some distance ^^o which is chosen to be much larger than 
the Vainshtein radius (typically ^cx) ~ 100). The boundary conditions in the DL case are given by the DL linear 
regime (fT5)) and read 



w{^ = 0) -- 


= 0, 


u{^ = 0) -- 


= 0, 


Wi^^ioo) -- 


2 

Soo 


«(^=^oo) = 


4 


3^00 


^(^=^oo) = 


= 0. 
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One may notice that the condition on w at infinity is imposed on w rather than on w. This is not necessary, but 
we found it easier to proceed that way in practice. Once the DL solution has been found, we can turn to the fully 
nonlinear system. In that case, the boundary conditions are given by the linear solution psp and read 

w{^ = 0) = 0, 
u{^ = 0) = 0, 

«^(C = Coo) = C X 1^ fl + ^ + ^^) e'^^-, 

3 4oo V «Coo «oo)V 

z^(C = Coo) = -Cx-- 



Pi^ = Coo) = 0. 

We have explicitly kept in these boundaries conditions the integration constant C defined in Eq. (|25p (while a is 
defined as in Eq. ([M]))- This procedure allows to check numerically that the choice we made C = 1 in Eq. (P5|) is in 
agreement with numerics. 

In our study of NLPF, we have implemented the relaxation method described in [2g|, using the C++ language. We 
typically used a grid of iV ^ 5 x 10"' points and aimed at a precision of 10^^. We always checked that our numerical 
solutions were stable under the change of the step size of the grid and the precision. 

The density of the star is assumed to be given by p = M/{4ttRq) inside the star, and null outside. 

B. Numerical results 

To summarize here shortly our main results, we were able to show that one can obtain a non singular solution 
showing the recovery of GR a la Vainshtcin. The relaxation method allowed us to solve the equations of motion ([5]) 
and ([9|) on a wide range of the parameter a = Rvm (a G [10~'^,0.6]) and for very large intervals (typically for R 
between and ^ 100 Ry)- In the following we first discuss the behaviour of the functions A and i' which appear in 
the "physical" metric g^i^ fsubsection I V B l]) . comparing their behaviour in NLPF and GR, and showing in particular 
how our numerical integration agrees with the corrections to GR as predicted by Vainshtein. Then (subsection IV B 2|) 
we turn to discuss the obtained behaviour of the pressure inside the star, comparing again the numerical solution 
with the one obtained in GR (recall that the star's density is assumed to be constant). Last we discuss our numerical 
result for the function ^ appearing in the "non physical" metric with our gauge choice ([7|). 

1. Behaviour of the functions X and v 
The Fig. [T] summarizes our main results and illustrates the following points: 



• 



at short distances {R <C i?v), the solution is very close to the GR solution, and in particular z/ ~ — A ^ —RsjR 
in the neighborhood of the source, 

at large distances (R ^ -Ry), the solution is in the linear regime (1^51) . Between Ry and m^^, the solution is 
very close to the DL solution (l?T|) . and in particular v ^ — 2A ^ —4/3 x RgjR. Beyond m~^, the solution is 
exponentially decreasing. 



• the behaviour at infinity allows to check that the choice of integration constant we made in Eq. p6| is consistent 
with numerics. More precisely, we tried to solve numerically the equations of motion with a integration constant 
C slightly different from 1, and we noted that the relaxation method was not converging, unless |C — 1| ^ 1. 

We have seen in subsection llll CI that the leading behaviour of the difference between the functions A and v and their 
GR equivalent, \gr and vqb, is given by the expressions (l47l) . It is possible to check that this theoretical prediction 
is very well satisfied by our numerical solution, as can be seen in Fig. [51 In this figure, we plotted the difference 
between the function A obtained from our numerical integration and the one of GR, and compared this with the first 
correction to GR as computed using the Vainshtein mechanism. As expected, the two agree very well for distances 
smaller than the Vainshtein radius. 

Besides, we have shown at the end of the subsection llll GI that deep inside the Vainshtein radius, the NLPF solution 
should be closer to the GR solution than to the DL solution, illustrating the fact that the GR nonlinearities are 
starting to play an important role. This is seen in Fig. [3J where we compared the function v (after an appropriately 
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FIG. 1: Plot of the metric functions —v and A vs. R/Rv , in the full nonlinear system and the decoupling limit (DL), with a 
star of radius Rq = lQ~^Rv and m x Rv = 10~^. For R <^ Rv, the rmmerical solution is close to the GR solution (where in 
particular f ~ —A ~ —Rs/R for R > Rq). For R ^ Rv, the solution enters a linear regime. Between Ry and m~^ , where the 
DL is still a good approximation, one has v ~ — 2A ~ —4/3 x Rs/R. At distances larger than m~^ the metric functions decay 
a la Yukawa as appearing more clearly in the insert. The latter shows the same solution but for larger values of R/ Rv, and in 
the range of distance plotted there, the numerical solutions are indistinguishable from the analytic solutions of the linearized 
field equations Eqs. (|25[l . This plot was already presented as the Fig.l of Ref. |26,] . 



chosen rescaling) for various values of the parameters. It appears on this figure that the numerical solutions well agree 
with the GR behaviour, while slightly differing form the DL expression for large a. This illustrate the fact that inside 
the source, the nonlinearities coming from the GR terms of Eq. ([8]), and which are not taken into account in the DL, 
are important. 

2. Behaviour of the pressure P 

In GR, the pressure inside a static spherically symmetric source of constant density can be computed analytically 
to be 



Pgr = p- 



l-^R2_ i_^ 



Rq 



3J1- 



Rs 

Rp, 



1 _ % i?2 



(61) 



This is compared on figure |3] with our numerical solution, confirming once again the validity of the Vainshtein's 
conjecture. 



3. The function fi 



For the range of parameters investigated here, the numerical behaviour of the function /i between R = and 



R 



is very close to the Vainshtein solution (|2ip of the DL. Beyond 



the Yukawa decay seen in the linear 



solution (1251) is found to occur. These various regimes are shown in Fig. [5] which plots the (rescaled) function /i for 
three different values of the parameters. Figure H] shows the same, zoomed inside the star. We notice there that the 
solution of the full theory goes towards DL solution as the parameter a = Ry x m tends towards and we recover 
the fact that the corrections to the DL due to the nonlinearities which are not taken into account in the DL are of 
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FIG. 2: Plot of the difTerence between the numerical solution for A, Xnum and the one of GR, Xgr- This is compared to the 

computation of the same quantity using the first term in the expansion in m^ proposed by Vainshtein, |AAy| — -^ (mR) \/ -^■ 

Both functions are plotted vs. R/Ry, for a source of radius Rq = 10~'^-Rv and a choice of parameters such that a = mx Rv ~ 
10~^. It can be seen that the numerical solution agrees very well with the Vainshtein's expression given by Eq. H47p for 
R <C Rv- For R ~ Rv, the expansion in m^ is, as expected, no longer a good approximation of the solution. 



order 0{a^) and goes to zero as a — ?■ 0. This plot also shows that our numerical integration reproduces nonlinearities 
not present in the DL. 



C. Shooting method 



In addition to the relaxation approach that we presented in the previous subsections, we developed a shooting 
approach to solve the system of equations ([5]) and ^ . This shooting approach, which is based on a direct Runge-Kutta 
integration, appears to be quite difficult to implement for the theory ([3]). Indeed, there appear numerical instabilities 
for distances beyond the Vainshtein radius and probably related to the growing mode ([55)1 . As a consequence of these 
instabilities, the initial conditions need to be extremely finely tuned^. Nevertheless, we have been able to integrate 
the system of equations outward from a point deep inside the star up to a point further than the Vainshtein radius 
(typically R ^ 3i?y). We have also been able to integrate inward starting from a point far from the source {e.g. 
R ~ 3i?T/) up to the core of the source. In all the integrations, we have made an extensive use of series expansions of 
the solution in order to fix the initial conditions as well as possible. 

While integrating inwards from large distances {i.e. R ^ Si?^), the solution adopts almost systematically a Q- 
scaling behaviour, similar to the one of Eq. (J19p . except if the initial conditions are extremely fine tuned to get 
the Vainshtein solution. In the cases in which the Q-scaling is found, our numerical solutions always encountered a 
singularity around the radius we identified in Eq. (j22p and which corresponds to the scale at which nonlinearities non 



"^ Note however that the shooting integration appears more stable for some NLPF theories which have a different interaction term from 
(|3]l. This is the case for example for the BD theory (using the terminology of Ref. 1291 ). In this theory, there is no admissible Vainshtein 
scaling, but only Q-scalings in the DL. However, in the DL and in the full non linear case, the shooting method appears much more 
stable than in the model considered here, even though the Q-scaling solution does not seem to continue into a non singular solution in 
the beyond DL regime. 
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FIG. 3: Plot of —1/ X a^* vs. R/Rv (the a~* factor is included for convenience such that, in the decoupling limit (DL), all 
plotted theories would exactly coincide) inside the source of radius Rq — W~'^Rv, for three different values of a = tti x Rv- 
Along with the numerical solution, the GR analytical solution for u is plotted for each value of a. It can be seen that the 
numerical solutions perfectly agree with the GR behaviour, while slightly differing form the DL expression for large a. This 
illustrate the fact that inside the source, the nonlinearities coming from the GR terms of Eq. (|8]), and which are not taken into 
account in the DL, are important. For small a, both the GR and the numerical solutions agree with the DL solution, which 
encodes for most of the physics, in agreement with the fact that the DL corresponds to the limit a — >■ 0. 



taken into account by the DL start to play a role. Of course, this does not mean that a Q-scaling solution cannot be 
found but we were not able to find any and it convinced us to focus on the Vainshtein solution. 

To obtain the latter, we have found particularly convenient to integrate from a point located inside the source 
(typically Rmin = IO^^Rq). To do so, we needed initial conditions at Rrmn, i-e. ^.{Rniin)-, fJ-'iRmin), K^min), 
v{Rmin) and P{Rmin)- Thcsc initial conditions are computed using the expansion (|49p where the three integration 

constants /.tg ' = /i(i? = 0), v^' = v{R = 0) and P^' = P{R — 0) have been fixed using the values obtained through 
relaxation method. In order to be able to obtain a wide enough range of integration, the required relative precision on 
the three integration constants /Xq 



W ,M 



from i?j, 



lO-^i?^ to R„ 



(x) 

and Pq is very stringent. For example, if one wants to be able to integrate 



3-Ry, in the case a = Rym = 10 ^ and Rq — 10 ^Ry, the required relative 



,(^) 



precision on the integration constant ^q is smaller than 10~ , in order to ensure that the numerical solution at 
Rmax = 3i?y differs from the linear solution (|25p (which is supposed to be an excellent approximation of the solution 
at that distance) by less than 5%. The result of such an integration is given in Fig. [T] This figure also shows part of 
the expansion (|^ inside the star. It agrees very well with the numerical solution. 



D. Weak- field approximation 



We have also been able to check numerically the weak field approximation introduced in section IIVI Namely, as 
appears on Fig. |S1 this approximation reproduces very well the numerical solution on a range of distances much larger 
than the decoupling limit, since this range extends at distances larger than the graviton Compton length m~^ . At 
small distances, however, as also seen on figure [51 the week field approximation deviates from the numerical solution. 
This is expected since the week field approximation reduces there to the DL and does not capture some non linearities. 
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FIG. 4: Plot of the functions P/p and Pgr/p vs. R/Rv, for a source of radius _R0 = 10~^Rv, of pressure P, and of constant 
density p, and for a parameters choice such that a = m x Rv ~ 10~^. The function P/p if found numerically, while Pgr/p is 
the pressure computed in GR and is given by the analytical expression (|6ip . The two curves are indistinguishable from each 
other, illustrating the fact that the Vainshtein's conjecture is valid inside the star. 



VI. CONCLUSIONS 



We have studied in detail the Vainshtein mechanism for Non Linear PauU-Fierz theory as defined in section|TIl More 
specifically, we have been able to show in the case of a specific NLPF model (the one defined by the interaction term 
([3])) that the recovery of the Schwarzschild solution works as forecasted by Vainshtein. This proof was provided via 
numerical integration of the field equation using two different methods: a relaxation scheme and a shooting method. 
It also heavily used series expansions (to fix boundary/initial conditions) which were first presented in this work. Our 
result, already presented in brief in [21| disagrees with a previous work on the same issue [201, which concluded that 
the Vainshtein mechanism was not working, because it could not find numerically any non singular asymptotically 
flat solution of the NLPF field equation. We also presented arguments, based on various expansions, showing that the 
solution at infinity is not uniquely defined by standard perturbation theory, but has non-perturbative hairs allowing to 
match sources at smaller distances. This, as well as the more sophisticated integration methods we used, could explain 
the differences between our results and the ones of Ref. [20| . Besides the model defined by the interaction term ([3]) , 
we also investigated other models, in particular those which were shown in Ref. [26] to possess only a non- Vainshtein 
like scaling (dubbed Q-scaling in Ref. [2g|) in the Decoupling Limit (DL). According to our numerical investigations, 
the DL solutions with Q-scaling, do not seem to continue into non singular solutions of the full field equations. We 
have also introduced a new limit, the weak-field approximation, which captures all the salient feature of the solution, 
including the Yukawa decay and the Vainshtein crossover. Of course, the theory considered is believed to suffer from 
various pathologies, like Boulware-Deser type of instabilities [3,|23,|30|- However, our investigations show for the first 
time that the Vainshtein mechanism can actually work in theories of massive gravity, some of which could be safe 
from the pathologies of NLPF theory. Several questions are however left unanswered. First, we have not been able 
to solve numerically the equation for dense objects (or even realistic stars). Indeed, when one increases the density 
of the object, the numerics becomes unstable and singularities are found to appear. There remains to understand 
if those singularities are physical or just numerical artifacts. Along the same line, it is not known if standard black 
holes can be recovered via the Vainshtein mechanism. Second, it is not clear if the solution we found is stable or not. 
The Boulware-Deser type of instability could show up if one would try to let the solution be time dependent or to be 
reached by a dynamical process such as spherical collapse. These and other issues are left for future investigations. 
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FIG. 5: Plot of a~^fi vs. R/Rv (the a~ factor is included for convenience such that, in the decoupling limit (DL), all plotted 
theories would exactly coincide) for a source of radius Rq = 10~'^Rv for three different values of a = m x Rv = 0.005, 0.05, 0.1 
as well as in the DL case (which corresponds to a — > 0). The three DL regimes of Eq. (|2H) are clearly distinguishable: for 
R < Rq, ^j. r^ const; for Rq < R <^ Rv, n ~ y^ {8Rv) / {9R); for Rv <, R <, m'^ , fi ~ 2i?^/(37?^). For R > m"\ we 
observe the Yukawa decay of the linear solution (|25p . At the resolution of the picture, the thee numerical solutions are hard 
to distinguish for distances R smaller than m~^. The next figure shows a zoom of this figure at small distances (namely inside 
the source) where the differences between the three solutions are easily seen 
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Appendix A: Details on perturbation theory 



Here, we explain with some details how one can obtain the series expansions (|28p and pip , and how the former is 
related to the DL expansion (|29| . Our starting point is the system ([8]) rewritten as in (|27| . 

Looking for an expansion of the form (|23p , and assuming such an expansion has been found up to order n— 1 (with 
71 > 1), the n th order system is of the form 



z 



z^ 2 



\i 



1 



z z^ 2 

^_ ^ - if 

~2 0-, ~ ^-/ff."' 



RR.,m 



(Al) 



z^ 2z z" 
where the functions /t^„, /ij_„, fg.n, Gtt,n and Gjm^n correspond respectively to the terms of order e"+^ of ft, Jr, fg, 
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FIG. 6: Plot of a~^fi vs. R/Rv (the o~^ factor is included for convenience such that, in the decoupling limit (DL), all plotted 
theories would exactly coincide) for a source of radius Rq — W~^Rv for three different values of a = m x Rv = 0.005, 0.05, 0.1 
as well as in the DL case. The picture is a zoom of Fig. [S] inside the star. One can check that even if the three curves for 
a = 0.005, 0.05, 0.1 are close to each other and to the DL solution, they differ inside the star due to non-linearities not captured 
by the DL. One note that the full solution goes towards the DL solution when a = Ry x m tends towards 0. 



Gtt and Ghr- It is possible to reorganize the system (|Aip to obtain the equivalent system 
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(A2) 



(A3) 



(A4) 



z^ z 

where L stands for the linear differential operator 

L[f] ^ -3z(4 + z')f{z) + 6{z^ + 4)/'(z) + 3ziz^ + 2)f"{z). 

Let's first focus on equation (jA4l) . which allows to find /^„. This equation can be written in the schematic form 

where F gathers the nonlinear terms of the equation of order e"+^, which only depend on the terms of order m < n. 
Note that the above discussion holds both for the large z and the small z expansions, since it is only based on the 
matching between terms of the same order in the parameter e on the left and right hand sides of equations (P7| . To 
go further, let us consider the action of the linear operator L on one term of the series (pS)) for /i„ of the form 
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FIG. 7: Plot of the function /i {fiN as given by the numerical integration) vs. R/Rv, for a source of radius Rq = 10 ^Rv and 
of constant density p, and for a choice of parameters such that a = m x Rv = lO"'^. The plot has been obtained using an 
direct integration d la Runge-Kutta. We also added the plot of series fj,s = /ig + fJ-i x'^ + fj.2 x'^, which was introduced in the 
main text as an expansion inside the star: the two curves are indistinguishable from each other. 



such that, with those notations, one is looking for an expansion of /i„ in the form 

oo n 
i=0 fe=0 

Henceforth, we will denote by Sn,m the set of formal series of the form 

n+l oo n 



(A6) 



Z'' 



^ ^ Sn,i,k z' {z^ log zY 



i=0 k=0 



where Sn,i.k denotes real coefficients. With this notation, we see that the expansion for /x„ given in Eq. (|28p and 
(|A6[) belongs to the set Sn,2, while the expansions for A„ and Vn belong to the set Sn,4- It is easy to check that those 
expansions hold for n = 0, as can be seen from equation (j25p . One finds easily that 

zL[Mn,i,k] = Pn,i,kMn^^^k 

+ Cln,i,kMn,i+l.k + 'r'n.i,kMn,i+2,k + Sn^i,kMn.i+s,^k + tn.i,kMn.i+4,k 



+ {k — n) [an.i,kMn.i+5,k+l + i>n.i^kMn,i+&,k+l + Cm,i,kMn^i+7 .k+1 + dnA^k^nA+S^k+l] 

+3(A: -n){k-n+ l)(2Af„,j+io,/c+2 + Mn,i+i2,k+2) 



(A7) 



where one has 



Pn.i,k = 6(i - 5fc)(z - 5fc 
0"n,i.k = 6(3 — 2i + lO/c) 



3) 



and qn,i,k, fn.i.k, Sn,i.k and tn.i,k are quadratic polynomials of n, i and k with integer coefficients, bn^i^k, Cn.i,k and dn.i,k 
are linear polynomial of n,i and k with integer coefficients. The exact expressions of these polynomials will not be 
needed here (except for tn.i,n, see below). It is then easy to see, that the action of the linear operator L on any series 
of Sn,2 results in a new series element of Sn,i- At each order n, a necessary condition to find a solution for /i„ in the 
above form (jA6p is then that F[{fim<n, ^m<m Vm<nY\ can be written as an element of Sn,i- Let us show that this is 
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FIG. 8: Plot of the function a^^fi (the a~^ is added to make the comparison between /i and its DL behaviour easier), vs. R/Rv, 
for a source of radius Rq = 10~^Rv, of constant density p and a choice of parameters such that a = mx Ry = W~^ . The three 
curves corresponds respectively to the weak- field solution of the equation (|59p . to the solution of the fully nonlinear theory and 
to the DL solution. One can see that the weak-field model allows to recover the general behaviour of the full solution, except 
inside the star (as seen in the panel showing a zoom of the plots for small R) where it is very close to the DL, and thus misses 
some nonlinearities. 



indeed the case. To do so, let us first consider the different functions appearing on the right hand side of Eq. (|27p . 
Those functions have different dependencies on A, v, /i and z that can be read from equations (|8all8bp and (llOalllOcj) 
and are given below. 

• the functions ft.^2 and fR^^2 come from the mass term and, as a consequence, can only depend on A, z^, /i, zfj,'. 
This means that once expanded in terms of powers of A, u, fj, and z^' , the functions /t.^2 ^ind fR.^2 are sums of 
/"2/i"^(z/i')"*, with ni,ri2 and n^, positive integers and n^ an integer verifying < n^ < 2. 



terms of the form A" 



• the functions Gtt.n and Grr^u come from the Einstein's tensor, in the gauge ([T]). Their expression is given by 
the left hand side of the equations ([5a|) and (j8b[) : one can then easily realize that the functions Gtt,n and Grr^u 
can be written as an infinite sum of terms of the form X""^ v"^ z~'^ , X"^i'"^i''z~'^ or X^^i'^^X'z~^, with rii and 
n2 positive integers. 



the function fg corresponds to the Bianchi identity 



to the derivative of the mass term; it is an infinite 



sum of terms of the form z ^A"^i^"^/x"^(z/i')"*(zz^')"'^(zA')"''(z^/i")"^ where Ui {i = 1...7) are positive integers 
and < 714 < 3, and < n^ < 1 for 5 < i < 7. 

It is then clear that, if one assumes that the series expansions hold at order e'' , with k < n, the right hand side of 
equation (jA4p can be written as a series of the form 






oo 



n 



k=k„ 



i — 'i-min k — k„ 



1 — fc 



m„,,,feM„,- fc = (6 e^)"-^ >J J2 m„,,,feZ^-3-5'=(logz)"-- (A8) 



It remains to show that fcmin = and imin — ^1 for the above series to be in the Sn,i set. It is clear that the first 
relation holds because the highest power of log z that can possibly appear on the right hand side of (IA4[) is necessarily 
strictly smaller than the power of e, i.e., n+1, this because it is so in each term of the series Sn,m- Let us now compute 
imin- To do SO, wc uccd to look at the terms on the right hand side of (|A4[) which are the most singular power of z 
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when z goes to zero. Keeping in mind that the series S'„_2 are more singular when z — >■ than the series S'„.4, looking 
at the above described decompositions of the functions /t,^2, /fl,>2, Gtt,ni Gnn^ni and fg, and inspecting Eq. (|A4p . 
one sees that the most singular power of z of the right hand sides of this equation, comes from the terms proportional 
to z~'^fg_n, z~^f'g,^, and /g'„; and more specifically from products of ^"*, (z^')"j and (z^^")"' appearing there (with 
ni,nj and ni positive integers strictly smaller than n). In those products, the piece containing the most singular 
power of z will be given by terms of the form* 



j=J 



-3 






.N. 



with Nj, fij, kj and J positives integers verifying kj < rij < n, J > 2 and 

j=J 

Y,inj + l)N,^n + l, (A9) 

where the above equality selects the order e"+^. The above expression (jA10[) can be rewritten as 

Defining then k and i as 

j=J 
k = ^(l + fcj)iVj -1 

, = -5 + 2^7V„ 

the above expression reads now 

(e e")"+i2»-3-5fc (log ^)"-fc _ (AlO) 

With the above definition, i takes its minimal value equal to —1 for ^2]=! ^j — 2, which is achieved for J = 2 and 
Ni = N2 — I (for which there is always a solution to equation (IA9|) '). This proves that the right hand side of Eq. 
(IA4[) is in the Sn,i set, as announced. Moreover, taking now into account the logarithm, it is easy to see what are 
the most singular terms in the right hand side of Eq. (jA4p (as z goes to zero). Those are obtained by the terms 
(|A10[) with i — —1 and k — n, which is achieved again under the same conditions as above: X]j=i Nj = 2, J = 2 and 
iVi = A^2 = 1 (with the additional condition that Vj, kj = rij). This shows in fact that those most singular terms are 
simply given by the serie expansion of the Decoupling Limit. Indeed, taking this limit, the only terms which are left 
in the right hand side of Eq. (IA4[) are those in fg which are quadratic in /^. The same reasoning also shows that the 
non linear pieces on the right hand side of Eqs. (jA2p and (jA3[) are respectively in the S'„.3 and 5„.4 families and have 
their most singular piece at small z given by the part of /g_„ which is quadratic in ^. 

This first shows that a (formal) solution of the system (|A2IIA4p can a priori be searched in the form of the series 
(P5)l . but also that, provided such a solution can be found, a subseries of it which corresponds to the most singular 
terms at small z will be given exactly by the series expansion obtained in the Decoupling Limit (Eq. (jlSI) ). 

Let us now show explicitly how to obtain the coefficients ^J.n,i,k, and solve Eq. (IA4[) for /i„. It is first of interest to 
look at the expression (|A7p for k — n which reads 

zL[Mn,i,n] = 6(i - 5n)(i - 5n - 3)M„^i,„ 

+ (6n + 3n^)Mn,^+4.u (All) 



The factor o{ z ''in front of the expression below conies from a factor ^ ^ in front of /g.„ and another factor of z ^ inside the expression 

of fg,„ 
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The necessity to incorporate logarithms in the series expansion for /x„ clearly appears above. Indeed, imagine we look 
for an expansion without logarithms, i.e. at each order n, we look for an expansion for /i„, Vn and A„ of the form 

oo 

1=0 

oo 

i=Q 

oo 

i=0 

Note that this form is equivalent to the form (|A6p for n = 0, and it is in agreement with the zero order solution given 
in equation (j25p . and we have proven above that those series must start from i = 0. Then, at fixed n, assuming the 
above form holds for terms of order & with k < n, and matching formally the terms proportional to z~^ x M„ ^ „ on 
the left hand side and right hand side of equation (jASp . one can try to solve for the values of the coefficients m„^i, 
starting from z = to increasing values of i. This is possible only until one reaches i = 5n where the coefficient 
in front of Mn^i^n in (|Alip vanishes. This vanishing means that the coefficient in front of the term z~^ x Mn,5n,n 
on the left hand side of (|A5P will not depend on /x„.5„. On the other hand, this coefficient will only depend on the 
coefficients /i„.i with i < 5n which are already determined by matching the terms Mn^i^n, with i < 5n, and one has 
thus no freedom left to match the term proportional to Af„^5„^„ on the two sides of (jASp . This means that it is not 
possible in general to find for /i„ a series expansion of the above form (jA12p . Note that the coefficient in front of 
Mn^i^n in (jAlip also vanishes for i = 5n + 3, and hence one has the same problem for finding the coefficient /in,5n+3- 
Including logarithms however solves this problem and allows to get a formal series solution in the family 5'„_2 as we 
know show. 

Indeed, let us fix n and see how one can compute the coefficients ^n,i,k- Let us first assume that all the coefficients 
Mn,i,fc, with k < K [K some chosen positive integer) and i < 5k are known, and see how to obtain the coefficients 
IJ^n,i,K, with i < 5K. This is immediate, indeed, the terms proportional to z~^ x Mn.i,K (for some i < 5K) on the left 
hand side of equation (jASP depend only on the coefficients /in,i,fc, with k < K and i < 5k, as can be seen from Eq. 
(|A7p . Hence, by matching formally those terms on the left hand side and right hand side of equation (jASP one obtains 
a system of 5K equations for the 5K unknown coefficients fJ-n,i,K, with i < 5K. This system is invertible because, as 
follows from equation (jA7p it is represented by a triangular matrix with diagonal elements, given by the Pn.i,K, which 
are non vanishing (this because we assumed i < 5K). Hence one can solve for the fin,i,K^ and by recursion one obtain 
all coefficients fJ.n,i,k^ with i < 5k. The coefficients fin,5k,k can then be determined as follows: fJ-n,o,o is determined 
uniquely by matching the term proportional to z~^ x Af„^5_i on the left hand side and right hand side of equation 
(|A5p . Indeed, /i„,5,i does not appear in this term, since Pn,5,i vanishes, but it is replaced there by /in,o,o because 
On, 0,0 does not vanish. By recursion, one can obtain in a similar way all coefficients fin,5k,k, with k < n, indeed, one 
can check that the coefficients an,i=5k,k — 18 and hence do not vanish, which means that at fixed K, fJ-n,5{K-i),K is 
uniquely determined by looking at the term z~^ x Mn,5K,K- The last coefficient /i„,5„.„ is left undetermined and can 
be chosen at will being related to a zero mode, as we will see below. The coefficients fJ,n,i,k, with 5fc < j < 5fc + 3 
are then determined in a similar manner as those with i < 5k. One then obtain the fin,i=5k+3,k in way similar to the 
way one obtained the fj.n,5k,k (checking again that an^i=5k+3,k = — 18 do not vanish). One is left with one arbitrary 
choice of fJ-n,5n+3,n, and the remaining coefficients iJ-n,i.k, with i > 5fc + 3 are obtained inverting an infinite triangular 
matrix (or a finite one, if one restricts oneself to some maximum value of i). The two arbitrary constants to be fixed 
each correspond to one of the two independent zero modes of the operator L. One of this zero modes is fiQ given in 
Eq. (I25p . another independent zero mode is given by 



1 1 \ , 1 / 1 1 
- P" - - 1 - - + -2 



f^QMs = -(l + - + ^)e "=--(1-- + ^ 



z \ z z^ z \ z z 



One can check that both functions /ip and ^q^ms can be expanded as in (|A6p with all coefficients ^n,i,k vanishing but 
those with k = n, and i > 5n for /iq and those with k — n, and i > 5n + 3 for /io.bjs- Note, as discussed in the 
main text, that the inclusions of such zero modes in the series expansions might be necessary to enforce the hierarchy 

Having shown how one can obtain /i„ in the announced form, it is easy to obtain the corresponding expansions for 
A„ and z^„ from equations (|A2p and (|A3p . There, there is no free constant left to be chosen. 
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Let us now turn to the asymptotic expansion (|3T|) for large z. Following a similar reasoning as previously, it is first 
easy to see that, provided such an expansion holds up to order e", then the right hand side of equation (lASp will be 
a series of the form 



But one also sees from equation (jAlip that this form is also the one resulting from the action of the linear operator L on 
the series for /x„ given in (PT|) . The expression (lAlip then shows that one can solve uniquely for all the coefficients [in.i 
starting from the largest values of z, and no free constants appear in this case because the coefficient in,i,n = 6n + iv? 
does not vanish. One can then easily obtain the expressions for i/„ and A„ in the form (|3ip with no extra free constant. 



Appendix B: Series inside the source 

The coefficients of the series expansions ([^^ can be found solving order by order in x^ ; the first orders read: 
/i(x) = [i-^^ + 



1 



/ (a;) (x) (x) o (x) (x) (x 



(x) , (x) 



2e"''o -e^o -5e^o -''o + 24e''« +^« 



7n 



2(3e^r..rfpr-4)~2e''r(p(^)_^ 






x^+0{x^) 



4 Va;?, 



(-1 + e^^^') (-e-^o^' + e^o^' + 2e'^o^'+^o^') m^) x^ + O(x^) 



K^) =^^"^ + 



1 -i/<"'-2u<"' 



(x) „ (x) 



- 2e2A'^"''m2 + 2e'^o^'+2Mo^' ("m^ + P(|^) + ^'j') ^^ + 0{x*) 



,(^) 



^0 



_e-^^'^2 ^ gM<^'^^2 _ 2^2^^' ^2 ^ 2e"^o^'+2'^o^' f m2 + P^^^ + A 



x2 + 0(a;4) 



Appendix C: Infinitely many solutions at infinity: a series expansion approach 

In this Appendix, we would like to complete the argument on the existence, at infinity, of infinitely many solutions 
of the system ([5]) having the same asymptotic behaviour given by the linear solution (I25p . Our starting point is section 
IIIIBl where we assumed that there exists a solution a = {A, v, p,}, such that 

(T — > (To for z — > oo. 

where ctq = {Aq, i'q, /io}, we looked for a solution a as small perturbations of the known solution a: 

(J = a + da, 

where 6a = {(5A, (5;^, 5/i}. We gave the leading behaviour of these solutions at Eq. (pij) . In the following, we would 
like to go beyond this leading behaviour analysis and find a consistent way of looking for the solution da through a 
series expansion. 
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1. Equations for 5a 

The first step is to write down the equations for Sa, that would generalize the leading order equations ((33 
Schematically, taking into account (j32[) . the system of equations ([5]) can be linearized around the solution a, 





dFt 
da 


a- da' 


BFr 

da 


/'^^ da' 


dfg 
da 


s da' 





5a' = 0, 
5a' = 0, 
5/i" = 0. 



(CI) 



It is important to note that in the linearized Bianchi identity, the last equation of (|Cip . we do not assume that the 
terms containing S^l" are negligible, so we keep them. We are interested in asymptotic solutions at infinity, thus we 
can substitute a^ instead of a in (jCip . since ct — )■ ctq at z — >■ oo. 

Further, since we are considering a range of distances where the solution is deep in the linear regime, we can safely 
keep only the lowest order in powers of ctq: while making an expansion of coefficients of equations (jCll) . This leads to 
the equations. 



5\' 5\ 1 ,t.^ ^c r ;\ 

h ^ + - ((5A + 3(5^ + z5^i') = be 



7T- -{5v + 25^i) = be- 



6z-^ z'^ 6z J \z^ 6z J 
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Y 1 1 2 1\ „ / 1 1 ^ 
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( ^ 2\ /I 1 2\ , 








+ l3; + 3J'^+l-^-z-3J^^ 








/^ 2 2\ /I 1 1\ , 








/ 1 


4)*"" 









(C2) 



One can of course check that keeping only the leading behaviour in the limit z — )■ oo leads to the equations (|33|) . 



2. Ansatz for 5ct 



The next step is to adopt an appropriate ansatz, which can be chosen of the form 

<5M = F^(z)exp(J-/(z)e^/2), 



<5A = -FA(z)yexp(J-/(z)e^/2 



(C3) 



5v = -F.(z)^yr7i exp (-^ - /(z) e-/2) 
where the function / is assumed to satisfy the equation 

3z3/2 y(^) 



/'(^) = 



2Y^(1 + ^) 



(C4) 
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The form of this ansatz wiU be justified later. The formal solution of Eq. (jC4p is given by 

l-z o +3/2 

f{z) = e--/2 / ^ e"^ dt. 

The exact choice of the lower integral bound zq does not matter so much, since it will just translate into an overall 
constant factor in Eqs. (jC3|) . An asymptotic expansion of the solution can be easily found to be 

/(.) = ^(3.-f + ^ + ...). (C5) 

Note that the leading behaviour of the ansatz (jC3p matches the asymptotic behaviour (|4ip . 

3. The equations for the F^'s 

Once the ansatz (jC3p has been chosen, we can compute the equations for the F^'s. They read^ 

Et ^ e[°^ + (Vbe-^/^) eI'^ + (6e-^) i?f ^ + {b^/^e^^^/^) E^f^ + {b^e-^^) e\^^ = 0, 
Er ^ Ef + (V6e-^/2) £;« + {be~^) Ef + (b^/\-^^/^) Ef = 0, (C6) 

Eb ^ 4°^ + (Vbe-^-'^) 4'^ + {be-^) 4^ + {b"'e-'^/^) Ef ^ 0, 



Note that we have multiplied the tt equation by a factor \/be ^' ■^ in order to have consistent form among the three equations. 
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with 



4°) [F. 

i?f ) [F, 
Ff ) [F, 

Ff ) [F, 

Fi°) [F. 

Ef m 

Ef [Fr 



Ef [F. 



E^b'' [F^ 



Ef [F. 



3z5/2 



4^/ITT 
1 



[Fa(z)-F^(z)], 



-[2{z^ + z + 12) Fa(z) + 4zFj;(z) - (12 + z)F^{z) - ^zF^z)] , 
1 



4yiv^ 
1 



= [-2z2(z + 3)Fa(z) + (z^ + 2z2 - 3z - 3) F^(z)] 



- 24^2 [(6^=* - 10z2 - 52z)Fa(z) - (Sz^ + 24z2)Fj;(z) 
+ (z3 + 2z2 + 21z - 3)F^(z) + (4z3 + gz^ - 12z - 12)F;;(z)] 
= ^^z-^ + zF,{z), 



^i[FA(z) + F,(z)-2F^(z)], 
.^^^,^i^_[(2z3 + 3z^ + z + 2)F.(z) 
-4(z2 + z)F;(z)-6(z3 + z2)F^(z)], 

[6z(l + z)Fa(z) + (1 + 5z - 8z2)F^(z) + 4(1 + z)F'^{z)\ 



I2z' 



9zV2 



, [(z-* + 5z3 + 6z2 + 4z + 2)F^(z)] 

Vz + 1 ' ^ 'J 



[2Fa(z) + F,(z)-3F^(z)] 



[-6(z2 + z2)Fa(z) + (z2 + z + 2) F,{z) ~ 4(z + z2)F^(z) 



24^z5(z + l) 
+ (36z + 30)F^(z) + 24(z + z^)Fl{z)\ , 

[2z {Az^ - 7z2 + 21z + 20) Fa (z) + (8z^ + 8z2)Fj; (z) 



48z'' 



+8z (2z2 + 3z + 3) Fj.(z) + (35z^ + 27z - 4)F^(z) 
+(8z2 - 24z - 16)F;(z) - 16(z2 - 16z)F;(z)] , 

[(lOz'' + 21z^ + 18z2 + 9z + 6) F^{z) 



(C7) 



-4z (2z3 + 5z^ + 6z + 3)f;(z)] . 

Note that we have used Eq. (|C4[) to eliminate /' and /"; surprisingly enough, the equations (jC7l) do not contain / 
either. 



4. Expansion in powers of Vb e ^' 

The next step is to solve for the lunctions F;. To do so, let's assume they arc of the form 

oo 

F,(z) = ^^(")(z) (Vbe-^/^' 



n=0 



with FJ^" (z) at most polynomial, i.e. 

i^(")(z) = 0(zP^"'), 
and solve the equations (IC2p order by order in v6e~^/^. 



(C8) 



(C9) 
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a. 0-th order A priori, out of the three equations (jC6p . one can get the three 0-th order equations for F^ 



(0) 



E, 



(0) 



F 



(0) 



= 0, 



where k — {t,R,B} and i — {A, i', /i}. However, only two of them are independent. To see this more algebraically, 
one can write the 0-th order equations in a matrix form 



where 



The vector ^Ke 



m.f(°) = 0, 



M: 




and F^") 



(0) 

A 

(0) 

y 

(0) 



2 



Gz' 



3/2 



belongs to the kernel of M^ , while the vectors 



/z+l 



Va = 



3zV2 





V, = 



4\/z+T 
3z=/2 

2 



can be used to get the equations 



V^.M.F 
Vl.M.F 



(0) 
(0) 



F(")(z)-^(°)(z) = 



(CIO) 



The fact that M is of rank 2 is a direct consequence of the ansatz (jC3|) , and in particular of the differential equation 
(jC4p satisfied by /. More precisely, if one does not assume anything on the function /, and then computes the 0-th 
order equations out of the ansatz (IC3|) . they are still of the homogeneous form N.F''^-' — 0, but with the matrix N 
given by 



/^^(/(z) + 2/'(z)) 







N 



-^z{fiz) + 2f{z))\ 



V 



Vby/ziz+l){fiz)+2f'{z)) 



_Sz 



-1 



Vb^/z{z+lj{f{z)+2f'{z)) b{z + l){f{z)+2f'{z)Y 



12z2 



1223 



/ 



In order to have N of rank 2 (otherwise F'"' is trivially 0), it is necessary to impose that its determinant is 0. If in 
addition we require / to be positive, the condition detN = translates into the differential equation (jC4p . 

To conclude, we have found that the 0-th order equations are not enough to fix the F^ s; however, they allow to 

X 

order to solve for i^^ '{z). 

b. First order The first order equations read 



understand the origin of the condition (|C4[) and enforce F^ '{z) — F^ \z) = F^ {z). We now have to go to next 



p(0)rp(l) 



(i)rp(o)i 



Er\Fl"]+E):>\F.^ 



0. 



(Cll) 



Let's first project this set of equations onto the Kernel vector ^ Ker, since we know this procedure will cancel the 
Ef\F'>'^\ contribution: 



V 



Ke 



E(i)[J^(°)] = 



sVfe 



2z(z + l) 



(6z + b)F'>^\z) + 4z(z + 1)K1°)'(z) 







(C12) 



where we have used the relations (jC10[) to replace F\, F^ by F^. The equation (jC12p can be straightforwardly solved, 
leading to 



Fr(z)^^(°)(z) = F(°)(z) = 



A 



z5/4(l + z)l/4' 



(C13) 
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with A being an arbitrary integration constant. A few remarks can be made at this point. First, the arbitrariness 
of the constant A exphcitly shows that out of a solution a of the system of equations ([5|), it is possible to find a 
family of infinite number of new solutions, with the same asymptotic behaviour at infinity. Second, we note that 
f/°' = O (z~^/^), in agreement with the condition ((C9l) . 

Let's turn now to the two other first order equations. Projecting Eq. (jClip onto Va and V^, leads to 



^(1) 



W -^//^W - ^/^i^r(^) - 0, 



" ^ ' '' ^ ' 6V^5(z + l) '^ ^ ' 6v/z5(z + l) ^ 



These equations allow to replace F| and f!) by i^/i in the equations of order 2. These order 2 equations will lead 

(1) (2) (2) (2) 

to an equation for i^/j , and two other equations relating F^ and F;) to F^ , and so on and so forth, 
c. n-th order The n-th order equations read 



By iteration, it is possible to prove that 



pin, 



(z) - Fj") (z) - F^") (z) = O (z-3/2+») . (C14) 



We can sketch briefly how the demonstration goes: after having solved the equation of order n, the functions F^ 
are assumed to be known, and F^\fI)^ are can be expressed as functions of {Fj^^ ,F^ }. Then, projecting the 
equations of order n + 1 onto ^Ker allows to find F/j" , leading to the scaling (jC14p . Projecting onto Va and y^, 
relates F^" and F^" to FfT , closing the iteration process. 

To summarize the findings of this section, we can say that we have found a class of solutions for the system of 
equations (|C2I) . These solutions are of the form (|C3p , where the functions Fj can be expanded into the series of Eq. 
(IC8[) . There are an infinite number of such solutions. This shows that out of a solution a of the system of equations 
([5]), it is possible to find an infinite family of new solutions, with the same asymptotic behaviour at infinity. This 
freedom is crucial to match the vacuum solution with the one with source. 

5. Numerical validation 

To conclude this study of the infinitely many solutions satisfying the same boundary conditions at infinity, we 
checked that our solution (jC3l) fits well the numerical solution of the system of equations (IC2[) . An example is shown 
in Fig. m 
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FIG. 9: Plot of the numerical solution for (5^((") for the choice of parameters b = 0.1 and 5^(C = 10) = 1. We 
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